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Summary 


This  MURI  research  has  investigated  the  problem  of  rotorcraft  brownout,  which  is  one  of  the  mil¬ 
itary’s  most  significant  problems  encountered  by  helicopters  and  tiltrotors  when  they  are  operated 
in  a  desert-like  environment.  Brownout  occurs  when  loose  sediment  material  sitting  on  the  ground 
is  uplifted  by  the  intense  rotor  wake  generated  by  the  rotor  system.  The  consequences  of  brownout 
is  a  degraded  visual  environment  (DVE)  for  the  pilot  as  well  as  the  possibility  of  motion  cue 
anomalies,  which  poses  a  serious  safety  of  flight  issue.  An  interdisciplinary  research  program  has 
been  undertaken  to  advance  the  understanding  of  the  fluid  dynamic  mechanisms  responsible  for 
brownout  and  DVE  conditions,  which  has  also  led  to  various  types  of  predictive  methodologies. 
While  classic  mechanisms  of  sediment  mobility  under  the  action  of  a  flow  over  a  surface  covered 
with  a  sediment  bed  have  been  shown  contribute  to  the  problem  of  brownout,  several  other  unique 
mechanisms  of  sediment  uplift  also  arise  because  of  the  unsteady,  three-dimensional  and  strongly 
vortical  nature  of  the  rotor  wake.  The  rotor  blade  tip  vortices,  in  particular,  have  been  shown  by 
means  of  laboratory  experiments  and  modeling  simulations  to  be  responsible  for  mobilizing  and 
uplifting  most  of  the  sediment  particles  from  the  ground.  Several  types  of  modeling  approaches 
have  been  developed  in  this  MURI  research,  ranging  from  those  based  on  inviscid  assumptions 
with  semi-empirical  models  representing  sediment  mobility  and  uplift,  to  sophisticated  two-phase 
computational  fluid  dynamic  methods  based  more  on  first  principles.  Such  computational  ap¬ 
proaches  have  also  provided  a  rigorous  basis  to  better  understand  the  possibilities  of  developing 
brownout  mitigation  strategies  from  the  perspective  of  both  rotor  design  and  flight  path  manage¬ 
ment.  This  final  report  summarizes  the  work  done  under  this  MURI  research  program,  which  has 
run  from  August  2008  to  July  2014. 
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Introduction 


This  comprehensive,  Multidisciplinary  University  Research  Initiative  (MURI)  program  has  ad¬ 
vanced  the  understanding  the  rotorcraft  brownout  problem.  Brownout  is  a  phenomenon  that  occurs 
during  rotorcraft  operations  near  the  ground  in  desert  environments,  and  is  characterized  by  a  heavy 
dust  suspension  that  forms  around  the  rotorcraft.  The  intensity  of  this  dust  cloud  is  often  to  the 
degree  that  it  can  obscure  the  pilots  vision.  Furthermore,  the  rapid  motion  of  the  dust  cloud  that  is 
usually  produced  in  brownout  conditions  can  lead  to  the  pilot  experiencing  motion  anomalies  such 
as  vection  illusions.  In  severe  cases  of  brownout,  the  pilot’s  can  experience  vertigo,  which  can  lead 
to  a  loss  of  control  of  the  rotorcraft.  Brownout  is  the  Air  Force  Special  Operations  Command’s 
most  significant  operational  problem,  with  more  than  half  of  all  human-factor  related  rotorcraft 
mishaps  during  landings  and  takeoffs  being  attributed  to  encounters  with  brownout  conditions. 

This  MURI  research  program  has  been  structured  to  provide  advances  in  the  physical  under¬ 
standing  of  rotorcraft  brownout,  and  has  also  been  focused  toward  the  development  of  a  compre¬ 
hensive  methodology  to  predict  and  understand  how  it  might  be  possible  to  mitigate  the  effects  of 
brownout.  While  the  problem  of  brownout  has  been  better  understood  by  performing  research  into 
the  underlying  fluid  dynamics  of  how  the  action  of  the  rotor  wake  near  the  ground  mobilizes  and 
uplifts  sediment  particles,  brownout  is  also  a  broader  problem  that  comes  under  the  classification 
of  flight  in  a  Degraded  Visual  Environment  (DVE).  The  research  goals  of  the  MURI  have  been 
achieved  by  better  understanding  rotor  and  airframe  aerodynamics  when  operating  in  ground  ef¬ 
fect,  with  particular  emphasis  on  understanding  the  detailed  dual-phase  nature  of  the  rotor-induced 


Figure  1:  Photograph  of  a  landing  helicopter  beginning  to  encounter  brownout  conditions. 
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particle  flow  near  the  ground.  In  particular,  details  of  the  fluid  dynamic  interactions  of  the  rotor 
wake  with  the  mobile  sediment  bed  have  been  researched  to  expose  the  mechanisms  responsible 
for  sediment  particle  uplift  as  well  as  the  resulting  transport  of  particles  around  the  rotor  and  the 
rotorcraft  as  a  whole. 

The  formation  of  the  brownout  dust  cloud  is  a  fundamentally  complex,  two-phase  fluid  dy¬ 
namics  phenomenon  with  its  origin  in  the  effects  of  the  energetic  rotor  wake  on  the  loose  sediment 
particles  lying  on  the  ground.  Factors  that  influence  the  mobilization  and  transport  of  the  resulting 
two-phase  flow  include  the  processes  of  particle  mobilization,  entrainment,  suspension,  bombard¬ 
ment  and  deposition.  If  the  dust  clouds  become  sufficiently  dense,  then  a  fluid  dynamic  coupling 
can  be  produced  to  the  degree  that  the  particles  can  modify  the  rotor  wake  flow.  However,  the 
understanding  of  the  brownout  problem  is  extremely  difficult,  in  part  because  the  dust  cloud  for¬ 
mation  is  affected  by  many  parameters  including  (but  not  limited  to):  rotor  disk  loading,  blade 
loading,  number  and  placement  of  rotors,  number  of  blades,  blade  twist,  blade  tip  shape,  fuselage 
shape,  as  well  as  sediment  type  and  surface  condition.  Some  types  of  rotorcraft  seem  to  be  par¬ 
ticularly  susceptible  to  brownout  problems  as  a  consequence  of  the  reingestion  of  suspended  dust 
through  the  rotor  disk  followed  by  bombardment  of  this  dust  back  onto  the  sediment  bed  and  the 
ejection  of  yet  more  dust. 

In  this  MURI  research  program,  detailed  measurements  of  the  fluid  dynamics  of  the  rotor  wake 
as  it  interacts  with  the  ground  and  with  a  mobile  sediment  bed  have  been  conducted  at  various 
geometric  scales,  and  in  both  air  and  water,  to  provide  fundamental  insight  into  the  boundary  layer 
flow  on  the  ground  and  the  mechanisms  of  particle  mobilization.  While  most  of  the  experiments 
have  been  performed  with  the  rotor  hovering  over  a  ground  plane  and  a  mobile  sediment  bed, 
some  more  recent  dual-phase  flow  experiments  have  been  conducted  with  a  rotor  translating  over 
the  bed  and  so  representing  an  actual  flight  maneuver.  In  this  latter  case,  the  bundling  of  rotor 
wake  vortices  on  the  ground  plane  ahead  of  the  translating  rotor  has  been  shown  to  be  a  primary 


Figure  2:  Schematic  showing  the  mechanisms  of  particle  mobilization  and  uplift  as  a  heli 
copter  encounters  brownout  conditions. 
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mechanism  that  affects  the  mobilization  and  uplift  of  sediment  particles  and  hence  the  formation 
of  brownout  dust  clouds. 

The  computational  tasks  in  the  MURI  effort  have  been  focused  on  more  accurately  capturing 
the  vortical  wake  system  generated  by  the  rotor  and  preserving  this  vorticity  as  it  is  convected 
towards  the  ground.  The  modeling  work  has  included  of  rotor  induced  particle  flows  as  well  as 
vortex-fuselage-tail  rotor  interacting  flows.  The  computational  work  in  this  MURI  has  also  been 
concerned  with  using  advanced  techniques  and  architectures  for  the  computational  simulation  of 
brownout.  For  example,  methods  have  been  developed  to  improve  the  computational  speed  of  the 
techniques  used  to  predict  particle  motion,  such  as  by  using  clustering  methods  and  fast-multipole 
methods,  and  also  with  the  use  of  heterogeneous  computing  architectures.  In  particular,  particle 
clustering  methods  have  been  explored  to  significantly  reduce  the  computational  time  in  the  sim¬ 
ulation  of  a  large  numbers  of  convecting  particles.  However,  it  has  been  shown  that  clustering 
methods  must  be  applied  judiciously  so  as  to  preserve  the  accuracy  of  the  resulting  dust  cloud 
simulations. 

The  validation  of  the  various  modeling  approaches  has  been  an  important  goal  of  this  MURI 
research  program.  While  much  validation  work  has  been  done  based  on  outcomes  from  labora¬ 
tory  measurements,  including  validation  of  the  mechanisms  of  vortex-induced  particle  motion,  the 
ultimate  challenge  is  to  predict  the  behavior  of  a  brownout  dust  cloud  generated  by  an  actual  he¬ 
licopter  during  a  landing  in  a  desert  environment.  Because  of  the  challenges  in  learning  about  the 
problem  and  developing  modeling  tools,  progress  toward  this  latter  end  has  taken  much  longer 
than  was  initially  anticipated,  mainly  because  it  needs  to  synthesize  understanding  and  modeling 
advancements  made  in  the  other  tasks,  including  pilot-in-the-loop  simulations  that  can  be  executed 
in  real-time.  The  recent  focus  has  been  on  developing  a  hybrid  approaches  that  take  advantage  of 
the  best  features  of  the  work  that  had  already  been  performed  under  other  tasks.  However,  a  lack 
of  measurements  of  the  brownout  phenomenon  at  full-scale  (i.e.,  with  actual  helicopters  flown  by 
pilots  in  actual  brownout  conditions)  continues  to  impede  the  validation  of  modeling  activities. 

Meantime,  the  the  MURI  work  has  progressed  toward  better  understanding  the  actual  brownout 
problem  by  using  various  simulations.  The  MURI  team  has  used  a  state  of  the  art  a  flight  simulator 
with  advanced  visuals,  which  has  been  used  to  render  and  display  brownout  dust  cloud  simula¬ 
tions  in  real  time.  A  physics-based  and  GPU-accelerated  brownout  simulation  methodology  has 
been  successfully  developed  to  represent  the  dynamic  motion  of  dust  particles  and  the  optical  at¬ 
tenuation  caused  by  dust  particles.  Two  new  methodologies  were  developed,  namely  a  dust  cloud 
illumination  methodology  to  calculate  the  optical  attenuation  caused  by  dust  clouds  that  affects  the 
visualization  of  dust  clouds  from  the  pilot’s  perspective,  and  a  methodology  to  render  individual 
dust  particles  on  computer  graphics.  The  calculations  of  velocities  induced  on  dust  particles,  their 
positions,  and  the  optical  attenuation  have  been  implemented  on  GPUs.  In  the  future,  this  system 
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may  help  to  verify  the  fidelity  of  the  dust  cloud  simulations  through  pilot  assessments,  and  also 
begin  to  explore  other  strategies  for  the  mitigation  of  brownout  through  flight  path  management. 

The  research  outcomes  in  this  MURI  have  shown  that  brownout  can  be  successfully  mitigated 
by  means  of  modifications  to  aspects  of  the  rotor  design,  but  this  may  not  be  a  practical  approach, 
i.e.,  in  practice  the  rotor  design  is  driven  and  constrained  by  many  other  factors  besides  brownout 
concerns.  The  research  work  has  also  shown  that  brownout  can  be  mitigated  by  appropriate  flight 
path  management  by  the  pilot,  i.e.,  where  pilots  fly  prescribed  airspeed  and  altitude  profiles  to 
minimize  the  severity  of  brownout.  However,  the  number  of  variables  involved  in  the  flight  path 
management  problem  mean  that  the  ideal  or  optimum  flight  paths  will  be  dependent  on  the  type 
of  rotorcraft  (e.g.,  single  main  rotor,  tandem,  tiltrotor,  etc.),  the  weight  and  center  of  gravity  of 
the  aircraft,  the  physical  nature  of  the  sediment  and  the  sediment  bed,  as  well  as  potentially  the 
different  flying  skills  of  each  pilot. 

Mitigation  may  also  help  to  reduce  another  problem  plaguing  rotorcraft  operating  in  desert 
environments,  which  is  the  rapid  abrasion  seen  on  the  rotor  blades.  Abrasion  leads  to  rotor  per¬ 
formance  degradation  as  well  as  high  maintenance  costs  and  significant  downtime  because  of  the 
frequent  replacement  of  rotor  blades.  With  the  ability  to  predict  with  growing  confidence  the  posi¬ 
tion  of  particles  with  respect  to  the  rotor  blades  and  airframe,  future  work  may  lead  to  models  that 
can  improve  the  design  of  more  abrasion  resistant  rotor  blades. 

This  final  report  summarizes  the  progress  made  in  the  MURI  effort,  primarily  over  the  last  two 
years,  and  also  indicates  profitable  avenues  for  future  research  into  the  brownout  problem.  Further 
information  on  the  research  efforts  can  be  obtained  directly  from  the  Principal  Investigator  or  from 
each  of  the  respective  task  leaders.  This  MURI  research  has  also  had  a  strong  educational  and 
training  component,  and  the  students  participating  in  the  program  have  been  identified  under  each 
task.  Eleven  students  have  successfully  completed  their  Ph.D.  degrees  under  the  support  of  the 
MURI  and  a  further  twenty  have  completed  their  M.S.  degrees;  several  students  have  completed 
both  their  M.S.  and  Ph.D.  degrees  under  this  MURI.  Most  of  these  students  are  now  pursuing 
further  research  in  government  labs  or  in  industry.  In  particular,  one  Ph.D.  student  from  the  pro¬ 
gram  (Joe  Milluzzo)  has  established  a  successful  career  as  a  professor  at  the  U.S.  Naval  Academy. 
Another  student  (Monica  Syal)  works  in  the  rotorcraft  industry  and  was  recently  publicized  by 
NVIDIA  (a  manufacturer  of  GPUs)  as  being  one  of  their  major  success  stories  in  use  of  GPUs 
for  scientific  simulation.  Several  other  students  from  the  MURI  program  are  now  working  at  U.S. 
Government  laboratories. 
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restructured,  or  are  funded  separately.  **  New  task. 


Task  Descriptions  and  Summary  of  Results 


The  following  pages  provide  summaries  of  the  research  in  each  of  the  MURI  tasks.  Further  infor¬ 
mation  can  be  obtained  from  the  Principal  Investigator  or  from  each  of  the  task  leaders. 

Part  1.  Fundamentals  of  Rotor  and  Airframe  Aerodynamics  in  Ground  Effect 
Operations: 

•  Task  1.1  -  Measurements  of  Rotor  and  Fuselage  Aerodynamics  in  Ground  Effect 

•  Task  1.2-  High  Fidelity  Vorticity  Generation  and  Preservation  in  Ground  Effect 

•  Task  1.3  -  Fuselage  Configuration  Effects  on  Rotorwash  and  Brownout 

Part  2.  Fundamentals  of  Particle  Suspension: 

•  Task  2.1  -  Non-Uniform,  Near-Bed  Flow  Field  with  Impinging  Rotorwash 

•  Task  2.2  -  Two-Phase  PIV  Near- Wall  Bounded  Turbulent  Flows 

•  Task  2.3  -  Fundamental  Two-Phase  Measurements  in  Brownout  Fluid  Mechanics 

•  Task  2.4  -  Large-Eddy  Simulation  of  the  Interaction  Between  Vortices  and  the  Ground 

•  Task  2.5  -  Two-Phase  Large  Eddy  Simulation  Using  the  Mesoscopic  Eulerian  Formalism 

•  Task  2.6:  Effects  of  Model  Scaling  on  Sediment  Transport 

•  Task  2.7:  Turbulent  Near-Bed  Flows  Under  Impinging  Vorticity 

Part  3.  Brownout  Synthesis  Mitigation  and  Validation: 

•  Task  3.1  -  Pilot-in-the-Loop  Brownout  Simulation  and  Mitigation  Strategies 

•  Task  3.2  -  Development  of  Efficient  Airborne  Sediment  Tracking  Algorithms 

•  Task  3.3  -  Understanding  Brownout  and  Developing  Mitigation  and  Control  Strategies 

•  Task  3.4  -  ABATE  Simulation  Framework  and  Validation 

•  Task  3.5  -  Computational  Considerations  in  the  Prediction  of  Brownout  Dust  Clouds 

•  Task  3.6  -  Blade  Tip  Actuation  for  Tip  Vortex  Modifications  for  Brownout  Mitigation 
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•  Task  1.2-  High  Fidelity  Vorticity  Generation  and  Preservation  in  Ground  Effect 
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Task  1.1 


Measurements  of  Rotor  and  Fuselage  Aerodynamics  in  Ground  Effect 

Investigator(s):  J.  Gordon  Leishman 

Institution/Department:  University  of  Maryland/Dept,  of  Aerospace  Engineering 
Graduate  Student(s):  Joseph  Milluzzo  Contact  email:  leishman@umd.edu 

Background  and  Technical  Challenges 

The  problem  of  brownout  is  caused  by  the  action  of  the  rotor  wake  on  the  ground  below  the 
rotor,  which  creates  unsteady  flows  and  turbulence  that  act  to  mobilize,  uplift  and  entrain  sediment 
particles.  In  particular,  as  the  blade  tip  vortices  trailed  from  the  rotor  interact  with  mobile  particles 
lying  over  the  ground,  they  cause  locally  high  shear  stress  on  the  bed  and  liberate  many  particles 
up  into  the  flow.  The  tip  vortices  also  generate  intense  regions  of  upwash,  which  uplift  the  particles 
well  away  from  the  bed  and  into  the  main  flow.  While  the  tip  vortices  are  the  dominant  features  of 
the  rotor  wake,  the  wake  is  also  comprised  of  wake  sheets  containing  vorticity  and  turbulence  that 
are  trailed  from  the  inboard  sections  of  the  blades.  These  sheets  can  introduce  significant  amounts 
of  turbulence  into  the  developing  flow  on  the  ground.  Ultimately,  a  detailed  understanding  of  both 
the  trailed  tip  vortices  and  wake  sheets  from  the  rotor  system  is  critical  to  understanding  the  flow 
at  the  ground  and  the  development  brownout  conditions. 

The  specific  goals  of  this  research  task  are  to:  1.  Measure  the  effects  of  rotor  and  airframe 
design  and  rotor  operational  conditions  on  the  blade  tip  vortex  structures,  including  the  turbulence 
and  Reynolds  stresses  in  the  flow  at  a  ground  plane  below  the  rotor;  2.  Characterize  and  under¬ 
stand  the  detailed  structure  and  evolution  of  the  wake  sheet  and  its  contributions  to  the  developing 
turbulent  flow  at  the  ground;  3.  Examine  the  effects  of  airframe  shapes  on  the  flow  below  the  rotor 
at  the  ground  with  a  view  to  understanding  how  an  airframe  may  participate  in  the  onset  of  sedi¬ 
ment  mobility.  The  outcomes  of  this  research  have  furthered  the  physical  understanding  of  rotor 
and  rotor/airframe  configurations  in  ground  effect  operation  and  have  provided  key  measurements 
needed  for  the  validation  of  computational  modeling  being  done  in  the  MURI  effort. 

Technical  Approach 

A  schematic  of  the  setup  used  for  the  rotor  experiments  conduced  in  this  task  is  shown  in  Fig.  1.1. 
Both  one-bladed  and  two-bladed  rotor  systems  have  been  used,  with  blades  of  radius  0.408  m 
(16  inches)  and  a  mean  chord  of  44.5  mm  (1.752  inches).  An  integrated  thrust  and  torque  balance 
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Figure  1.1:  Schematic  of  the  setup  showing  the  rotor,  body,  and  the  ground  plane. 

is  used  to  make  rotor  performance  measurements.  The  rotor  system  has  a  variable  collective  pitch 
capability  so  that  the  measurements  can  be  compared  at  the  same  (or  equivalent)  rotor  operating 
conditions. 

In  the  first  phase  of  the  work,  the  effects  of  blade  tip  shape  on  the  flow  was  examined.  In  recent 
work,  the  effects  of  blade  twist  on  the  development  of  the  rotor  wake  was  studied.  Two  blade 
sets  were  tested:  1.  A  baseline  untwisted  blade,  and  2.  A  twisted  blade  with  —17°  of  linear  twist 
from  root  to  tip.  Each  rotor  system  was  tested  at  two  blade  loading  coefficients,  Cj / o,  of  0.053  and 
0.08,  which  were  selected  to  represent  helicopter  rotor  operations  at  low  and  high  thrust  conditions, 
respectively.  Flow  measurements  have  been  made  near  the  rotor  and  also  near  the  ground  plane. 

To  simulate  ground  effect  conditions,  a  circular  plane  of  diameter  equal  to  two  rotor  diameters 
was  placed  parallel  to  the  rotor  plane,  as  also  shown  in  Fig.  1.1.  The  rotor  was  located  one  radius 
above  this  ground  plane.  Flow  diverters  were  placed  around  the  periphery  of  the  ground  plane, 
which  controlled  the  flow  as  it  exited  off  the  surface  and  helped  to  limit  flow  recirculation  in  the 
test  cell.  Honeycomb  flow  conditioning  screens  were  located  upstream  to  reduce  turbulence  levels 
into  the  rotor  disk. 

Three  body  shapes  have  also  been  tested  in  the  wake  below  the  rotor;  see  Fig.  1.2.  These 
shapes  are  geometrically  defined  in  terms  of  super-ellipsoidal  functions.  The  bodies  were  mounted 
such  that  their  vertical  and  longitudinal  positions  with  respect  to  the  rotor  can  be  adjusted.  A  force 
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Figure  1.2:  Side  view  of  the  three  bodies  used  in  the  present  experiments.  All  shapes  defined 
in  terms  of  super-ellipsoidal  functions. 

balance  was  used  to  measure  the  download  on  the  bodies  under  the  action  of  the  rotor  wake. 

Various  measurements  have  been  conducted  with  the  rotor  and  each  body  shape  in  turn,  in¬ 
cluding  flow  visualization  (FV)  and  particle  image  velocimetry  (PIV),  as  well  as  performance  (i.e., 
thrust,  power,  and  body  download)  sweeps.  The  same  combination  of  lasers  and  camera  systems 
were  used  for  the  FV  and  PIV  experiments.  The  laser  was  fired  onto  a  high-intensity  rated  mirror 
and  then  into  spherical  and  convex  lenses,  which  produced  a  thin  light  sheet  in  the  regions  of  inter¬ 
est  (ROI).  The  camera  was  oriented  so  that  the  imaging  axis  was  orthogonal  to  the  laser  sheet,  and 
was  focused  to  image  the  desired  locations  in  the  rotor  wake  and  at  the  ground  plane.  The  camera 
and  laser  was  digitally  synchronized  such  that  the  laser  pulses  straddled  the  camera  images.  For  the 
FV  and  PIV  measurements,  the  camera  and  laser  were  synchronized  with  a  phase-locking  encoder. 

Results 

Figure  1.3  shows  the  regions  of  interest  (ROIs)  where  PIV  was  performed  for  the  isolated  rotor 
experiments  with  the  rotor  operating  out  of  ground  effect  and  in  ground  effect.  Multiple  interroga¬ 
tion  regions  were  needed  to  provide  the  necessary  spatial  resolution.  For  the  rotor  operating  out  of 
ground  effect  (i.e.,  ROIs  1  and  2)  and  in  ground  effect  (i.e.,  ROIs  1-6),  the  cameras  were  placed 
next  to  each  other  with  a  60  mm  overlap  in  their  field  of  view,  which  allowed  the  temporally 
correlated  images  to  be  mosaiced  together.  Each  camera  had  an  initial  ROI  of  270-by-180  mm 
(0.66-by-0.44  R ),  yielding  a  total  mosaiced  field  of  view  of  480-by-180  mm  (1.18-by-0.44  R ).  For 
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Figure  1.3:  Regions  of  interest  (ROI)  used  for  the  flow  measurements:  (a)  Out  of  ground 
effect;  and  (b)  In  ground  effect. 

the  rotor  operating  out  of  ground  effect  a  single  camera  was  focused  on  a  smaller  region  (i.e., 
ROI  3)  at  the  blade  tip  (approximately  102-by-68  mm  or  0.26-by-0.17  R),  which  provided  high 
resolution  images  (>  3,000  measurement  points  per  unit  radius)  of  the  developing  tip  vortices  and 
vortex  sheets. 

Figures  1 .4  and  1 .5  shows  representative  FV  images  that  were  taken  with  the  reference  blade 
at  yb  30°.  The  tip  vortices  produced  by  the  untwisted  blades  rolled-up  almost  immediately 
behind  the  trailing  edge  of  the  blade  tip,  with  the  vortices  effectively  being  fully-formed  at  wake 
ages  as  early  as  C,  =  3°.  However,  the  vortices  produced  by  the  twisted  blades  were  seen  to  take 
somewhat  longer  to  roll-up,  i.e.,  not  until  C,  ~  15°-20°  could  they  considered  as  fully  developed. 
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(a)  Untwisted  blade,  Ct/g  =  0.053  (b)  Untwisted  blade,  Cj / o  —  0.080 


Figure  1.4:  Flow  visualization  in  the  rotor  wake  near  the  tip  of  the  untwisted  blades  for  each 
of  the  operating  thrusts  at  a  blade  angle  of  30°:  (a)  Cj /a  =  0.053;  (b)  Cj /a  =  0.080. 


(a)  Twisted  blade,  Ct/o  =  0.053  (b)  Twisted  blade,  Cj/o  =  0.080 

Figure  1.5:  Flow  visualization  in  the  rotor  wake  near  the  tip  of  the  twisted  blades  for  each  of 
the  operating  thrusts  at  a  blade  angle  of  30°:  (a)  Cr/o  =  0.053;  (b)  Cr/o  =  0.080. 


Such  effects,  at  least  partly,  reflect  the  different  loading  distributions  over  the  tip  region  of  each 
blade.  Furthermore,  it  is  clear  from  Figs.  1.4  and  1.5  that  the  vortex  sheets  produced  by  the  twisted 
blades  were  significantly  thicker  and  more  turbulent. 

Figures  1 .6  and  1 .7  shows  the  measured  locations  of  the  tip  vortices  and  the  vortex  sheets  as 
a  function  of  wake  age  for  each  rotor.  A  wave  was  seen  to  form  on  the  sheet  not  long  after  it 
was  initially  trailed  from  the  blade.  This  wave  was  more  pronounced  in  the  wake  produced  by 
the  twisted  blades,  as  shown  in  Figs.  7(a)  and  7(b).  Moving  inboard,  this  wave  can  be  seen  to 
initially  grow  in  magnitude.  However,  as  the  vortex  sheet  convected  further  axially  in  the  wake 
it  was  stretched,  which  appeared  to  suppress  any  further  growth  in  wave  amplitude.  Additionally, 
the  sheet  trailed  from  the  rotor  with  untwisted  blades  was  seen  to  take  on  two  distinct  slopes  (i.e., 
negative  inboard  and  positive  outboard)  as  shown  in  Fig.  1.4;  this  effect  occurred  around  0.80/? 
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(a)  Untwisted  blade,  Ct/o  =  0.053 


(b)  Untwisted  blade,  Ct/o  =  0.080 


Figure  1.6:  Axial  and  radial  locations  of  the  wake  sheets  and  tip  vortex  locations  along  the 
wake  boundary  for  the  rotor  with  untwisted  blades:  (a)  Cr/o  =  0.053;  (b)  C7/0  =  0.080. 


(a)  Twisted  blade,  Cj/o  —  0.053  (b)  Twisted  blade,  Ct/g  =  0.080 


Figure  1.7:  Axial  and  radial  locations  of  the  wake  sheets  and  tip  vortex  locations  along  the 
wake  boundary  for  the  rotor  with  twisted  blades:  (a)  Cr/o  =  0.053;  (b)  C7 / o  =  0.080. 


at  each  of  the  two  rotor  thrust  conditions;  see  Fig.  1.6.  Notice  that  at  low  rotor  thrust,  the  tip 
vortices  produced  by  the  twisted  blades  were  convected  axially  more  slowly  than  for  the  other 
cases.  This  slower  convection  caused  the  following  blade  to  interact  with  the  vortex.  Such  head- 
on  interactions,  which  can  be  classified  as  form  of  perpendicular  blade  vortex  interaction  or  BVI, 
caused  the  rapid  diffusion  of  vorticity,  allowing  the  vortex  to  persist  in  the  flow  for  less  than  one 
rotor  revolution. 

The  locations  of  the  tip  vortices  and  wake  sheets  as  a  function  of  wake  age  are  shown  in 
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Fig.  1.8  for  the  rotor  with  the  untwisted  blades  when  it  was  operating  both  out  of  ground  effect 
and  in  ground  effect.  The  results  showed  that  the  non-uniform  convection  of  the  wake  sheets  was 
qualitatively  similar  for  both  operating  conditions,  with  the  quickest  axial  convection  occurring 
in  the  wake  outboard  of  r/R  =  0.6;  see  Fig.  1.8.  However,  the  behavior  of  the  wake  sheets  were 
significantly  affected  by  the  presence  of  the  ground.  After  C,  =  120°,  the  wake  sheets  from  the 
untwisted  blades  convected  axially  slower  below  the  rotor  when  it  was  operating  in  ground  effect. 
For  example,  Fig.  8(a)  shows  that  at  a  wake  age  of  C,  =  180°  only  the  inboard  75%  of  the  wake 
sheet  was  above  z/R  =  —0.35  for  the  rotor  operating  out  of  ground  effect.  However,  when  the  rotor 
was  operated  in  ground  effect,  the  entire  wake  sheet  remained  above  z/R  =  —  0.35  until  a  wake  age 
of  C,  >  210°  was  reached;  see  Fig.  8(b).  In  this  case,  the  slower  axial  convection  was  because  of 
the  presence  of  the  ground,  which  modified  the  spanwise  downwash  distribution  in  the  wake  below 
the  rotor. 

When  the  rotor  was  operated  out  of  ground  effect,  a  wave-like  perturbation  formed  on  the  wake 
sheets  shortly  after  they  were  trailed  from  the  untwisted  blade;  see  Fig.  1.6.  However,  as  the  sheet 
convected  axially  below  the  rotor  it  was  stretched,  which  effectively  suppressed  any  additional 
growth  in  the  wave  amplitude.  From  an  examination  of  Fig.  8(b),  it  can  be  see  that  when  the  rotor 
was  operated  in  ground  effect,  no  wave-like  perturbation  was  observed  to  form  on  the  sheet.  In 
this  case,  the  radial  stretching  of  the  wake  because  of  the  presence  of  the  ground  appeared  to  be 
the  source  of  the  change  in  the  behavior  of  the  wake  sheet. 

Figure  1.9  shows  the  locations  of  the  tip  vortices  and  wake  sheets  as  a  function  of  wake  age 
for  the  rotor  with  twisted  blades  when  it  was  operating  both  out  of  ground  effect  and  in  ground 
effect,  respectively.  Similar  to  that  of  the  rotor  the  with  untwisted  blades,  the  wake  sheets  trailed 
from  the  rotor  with  twisted  blades  convected  axially  slower  than  the  corresponding  out  of  ground 
effect  case.  Furthermore,  the  presence  of  the  ground  also  resulted  in  the  suppression  of  the  wave¬ 
like  perturbation  on  the  wake  sheets  trailed  from  the  rotor  with  the  twisted  blades,  as  shown  in 
Fig.  9(b). 

Imaging  of  the  tip  region  allowed  the  developing  wake  sheet  to  be  examined  in  detail.  Specif¬ 
ically,  the  near- tip  region  contained  over  70,000  PIV  measurement  points,  corresponding  to  over 
3,000  measurement  points  per  unit  radius.  Figure  1.10  shows  an  example  of  three  instantaneous 
contours  of  vorticity  (with  the  corresponding  velocity  vectors  being  superimposed)  at  the  higher 
thrust  condition  and  a  wake  age  of  6°  (i.e.,  C,  =  6°).  These  figures  show  that  the  wake  sheets  were 
comprised  of  concentrated  regions  of  counterrotating  vorticity  (often  called  Taylor-Gortler  vor¬ 
tices)  that  generated  significant  fluctuations  in  the  total  velocity  in  the  wake.  However,  the  twisted 
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(a)  Out  of  ground  effect,  (z/R  =  °°)  (b)  In  ground  effect,  z/R=  1.0 

Figure  1.8:  Plot  of  the  axial  and  radial  locations  (using  a  rotor  based  reference  frame)  of  the 
wake  sheets  and  tip  vortices  in  ROIs  1  and  2  using  a  rotor  based  reference  frame  for  the  rotor 
with  untwisted  blades  operated:  (a)  Out  of  ground  effect  (z/R  =  °°);  (b)  In  ground  effect  (z/R 
=  1.0). 
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(a)  Out  of  ground  effect,  ( z/R  =  00 )  (b)  In  ground  effect,  z/R  =  1.0 


Figure  1.9:  Plot  of  the  axial  and  radial  locations  (for  a  rotor  based  reference  frame)  of  the 
wake  sheets  and  tip  vortices  in  ROIs  1  and  2  for  the  rotor  with  the  twisted  blades  operated: 
(a)  Out  of  ground  effect  (z/R  =  °°);  (b)  In  ground  effect  (z/R  =  1.0). 


28 


Normalized  radial  distance,  r/  R 


Normalized  radial  distance,  r/  R 


(a)  Instantaneous  realization  one 


(b)  Instantaneous  realization  two 


Normalized  radial  distance,  r/  R 
(c)  Instantaneous  realization  three 

Figure  1.10:  Contours  of  instantaneous  vorticity  at  C,  =  6°  for  the  rotor  with  twisted  blades 
operating  at  the  higher  thrust  condition. 

blades  produced  a  more  intense  wake  sheet  that  contained  greater  vorticity,  which  also  leads  to 
more  turbulence  in  the  sheet.  The  turbulence  contributions  from  the  wake  sheets  are  of  particular 
interest  because  Task  2.7  has  shown  that  particle  mobilization  and  uplift  can  be  correlated  to  the 
Reynolds  stresses  and  discrete  turbulence  events  at  the  ground.  While  the  wake  sheets  appear  to 
contribute  significant  turbulence  into  the  developing  rotor  wake,  it  was  very  difficult  to  actually 
quantify  the  turbulence  because  of  the  highly  aperiodic  nature  of  the  wake  sheet.  These  instan¬ 
taneous  vorticity  realizations  show  the  inherently  aperiodic  nature  of  the  wake  sheet,  with  each 
image  showing  slightly  different  amounts  and  distributions  of  vorticity.  Because  the  size,  shape, 
and  structure  of  the  wake  sheet  varied  between  phase-resolved  realizations  at  the  same  wake  age, 
a  turbulence  analysis  cannot  be  performed  without  collocating  the  sheet  structures,  otherwise  the 
calculated  turbulence  values  would  be  biased  by  the  effects  of  aperiodicity. 

While  an  examination  of  the  flow  field  generated  near  the  rotor  provided  significant  insight  into 
the  development  and  initial  formation  of  the  flow  features,  the  problem  of  is  driven  by  the  near-wall 
flow  field.  Specifically,  the  problem  of  brownout  is  linked  to  the  structure  of  the  near- wall  vorticity. 
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Figures  1.11  and  1.12  show  representative  instantaneous  realizations  of  the  vorticity  in  ROIs  5  and 
6  for  both  blade  sets  at  a  reference  blade  azimuth  angle  of  \\ft,  =  30°.  A  considerable  amount  of 
small  scale  vorticity  in  the  wakes  was  observed  inboard  near  the  rotor  (i.e.,  for  r/R  <  0.6),  and  it 
would  be  expected  that  the  small  scale  vorticity  seen  in  Figs.  1.11  and  1.12  was  conveyed  there  by 
the  sheets.  Furthermore,  because  the  twisted  blade  initially  trailed  a  more  intense  wake  sheet,  the 
vorticity  in  the  near-wall  flow  would  be  expected  to  have  higher  magnitudes  as  compared  to  the 
corresponding  untwisted  blade,  as  shown  in  Figs.  1.11  and  1.12. 

As  discussed  further  in  Task  2.1,  the  formation  of  a  brownout  cloud  results  from  the  interaction 
of  the  tip  vortices  with  the  saltation  layer  along  the  ground.  Specifically,  the  mobilized  sediment 
along  the  ground  becomes  “trapped”  in  the  higher  velocities  induced  by  the  vortices  and  is  sus¬ 
pended  around  the  rotor.  Therefore,  the  question  becomes,  how  does  the  vorticity  in  the  near-wall 
flow  affect  the  wall-normal  velocities? 

To  help  answer  this  question,  contours  of  the  wall-normal  velocity  in  ROIs  5  and  6  are  shown  in 
Figs.  1.13  and  1 . 14  for  the  untwisted  and  twisted  blades,  respectively.  The  tip  vortices  are  identified 
in  Figs.  1.13  and  1.14  by  the  red  (upwash)  and  blue  (downwash)  pairs.  While  the  tip  vortices 
induced  strong  upwash  velocities  and,  therefore,  would  have  the  potential  to  uplift  large  amounts 
of  previously  mobilized  sediment  particles,  the  smaller  scale  vorticity  (shown  in  Figs.  1.11  and 
1.12)  also  induced  significant  upwash  velocities,  and  closer  to  the  rotor;  see  Figs.  1.13  and  1.14. 
It  would  be  expected  that  the  fluctuating  velocities  associated  with  the  wake  sheets  could  uplift 
sediment  particles,  and  so  will  be  a  contributing  mechanism  to  the  problem  of  brownout.  Without 
the  flow  fluctuations  induced  by  the  wake  sheet,  it  is  likely  that  the  saltation  layer  would  develop 
further  downstream  and,  therefore,  the  initial  uplift  of  particles  would  be  further  away  from  the 
rotor. 

A  body  representing  an  airframe  was  introduced  into  the  rotor  wake  flow  to  determine  its  effect 
on  the  developing  flow  over  the  ground  plane.  Scaled  down  bodies  have  also  been  tested  in  the 
dual-phase  environment  in  the  dust  chamber;  see  Task  2.2.  The  circular  and  elliptical  body  shapes 
were  found  to  produce  similar  flows  near  the  ground.  These  bodies  also  showed  near- wall  velocity 
and  turbulence  profiles  similar  to  those  obtained  with  the  isolated  rotor,  although  some  differences 
did  appear  at  the  most  inboard  radial  locations  nearer  to  the  rotor;  in  particular,  the  rectangular 
body  produced  different  Reynolds  stresses  at  locations  inboard  of  1.3 R.  These  differences  were 
particularly  significant  at  about  one  rotor  radius  where  the  wall-jet  had  not  yet  fully  developed. 
The  results  showed  that  at  downstream  distances  greater  than  about  r/R  =  1.6  there  were  minimal 
differences  in  the  flows  at  the  ground  as  influenced  by  any  of  the  bodies. 
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Figure  1.11:  Instantaneous  realization  of  the  flow  field  in  ROIs  5  and  6  for  the  rotor  with  the 
untwisted  blades  at  a  blade  azimuth  angle  \\ft,  =  30°.  Background  contours  are  of  vorticity. 
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Figure  1.12:  Instantaneous  realization  of  the  flow  field  in  ROIs  5  and  6  for  the  rotor  with  the 
twisted  blades  at  a  blade  azimuth  angle  \| ft,  =  30°.  Background  contours  are  of  vorticity. 


£  0  45 

N  04 


03  - 


25 

-Q 

m  0  2S 


0  2  ~ 


I  01$ 

I  °'f 

O  €05 


m 


n 


V.  ^  as  4  -15  D  ts  l  as  a 


03 


0.4  Ofi  Ofl  1 

Normalized  radial  distance,  r  f  R 


1.3 


t  4 


Figure  1.13:  Instantaneous  realization  of  the  flow  field  in  ROIs  5  and  6  for  the  rotor  with 
the  untwisted  blades  at  a  blade  azimuth  angle  \|/^  =  30°.  Background  contours  are  of  the 
wall-normal  velocity. 
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Figure  1.14:  Instantaneous  realization  of  the  flow  field  in  ROIs  5  and  6  for  the  rotor  with  the 
twisted  blades  at  a  blade  azimuth  angle  \\ft>  30°.  Background  contours  are  of  the  wall-normal 
velocity. 
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Task  1.2 


High  Fidelity  Vorticity  Generation  and  Preservation  in  Ground  Effect 

Investigator(s):  James  Baeder 

Institution/Department:  University  of  Maryland/Dept,  of  Aerospace  Engineering 
Graduate  Student(s):  Tarandeep  Kalra 
Contact  email:  baeder@umd.edu 

Background  and  Technical  Challenges 

The  action  of  the  rotor  wake  on  loose  sediment  on  the  ground  is  primarily  responsible  for  inducing 
the  rotorcraft  brownout  phenomenon.  Therefore,  any  simulation  of  brownout  must  be  capable  of 
accurately  predicting  the  velocity  field  induced  by  the  rotor  when  it  is  operating  in  ground  effect. 
This  goal  requires  the  computational  methodology  to  be  able  to  capture  the  generation  of  vorticity 
on  the  blade  surfaces,  including  the  tip  vortex  and  the  formation  of  the  wake  sheets,  as  well  as  the 
subsequent  evolution  and  rollup  into  the  tip  vortices,  and  the  convection  and  interaction  of  of  this 
vorticity  into  the  developing  wall  jet  and  boundary  layer  on  the  ground.  However,  the  accurate 
prediction  of  the  wake  vorticity  from  the  rotor  plane  to  the  ground  plane  with  the  use  of  high- 
fidelity  CFD  methodologies  is  computationally  expensive,  mainly  because  of  the  fine  grid  sizes 
that  are  needed  to  prevent  the  non-physical  numerical  dissipation  of  the  vorticity  at  older  wake 
ages.  The  increase  in  the  physical  scale  of  the  rotor  system  (i.e.,  from  laboratory  or  model-scale  to 
full-scale)  further  increases  the  computational  challenges  because  of  the  larger  sizes  of  the  grids 
that  are  needed  to  solve  for  the  flow  properties,  and  also  the  computational  time  needed  in  this  case 
can  become  increasingly  prohibitive. 

Technical  Approach 

Experimental  studies  conducted  in  this  MURI  (see  Tasks  1.1  and  2.2)  have  performed  detailed 
measurements  of  the  flow  near  and  below  sub-scale  rotors  operating  over  a  ground  plane.  These 
measurements  have  provided  valuable  data  to  help  validate  the  CFD  models.  After  validation,  the 
CFD  models  can  then  be  used  to  provide  results  that  can  expose  the  relative  sensitivities  of  various 
rotor  (e.g.,  the  blade  tip  shape)  and  particle  parameters  (e.g.,  the  particle  diameter,  mass,  etc.)  to 
the  severity  of  the  brownout  phenomenon.  Such  CFD  models  can  also  simulate  combinations  of 
rotor  and  particle  parameters  that  may  not  easily  be  studied  experimentally  in  the  laboratory. 

In  this  task,  the  overset  version  of  Transonic  Unsteady  Rotor  Navier-Stokes  (UMTURNS) 
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(a)  Micro-rotor  (b)  Sub-scale  rotor 

Figure  1.15:  Comparison  of  CFD  (S-A  model)  predicted  eddy  viscosity  contours  at  an  az 
imuthal  plane  of  0°.  (Notice  the  different  contour  levels). 


is  being  used  to  help  computationally  analyze  the  rotorcraft  brownout  problem.  A  prerequisite, 
however,  is  to  predict  the  flow  without  the  particles.  In  this  solver,  the  inviscid  terms  are  discretized 
using  a  third  order  MUSCL  scheme  with  Roe’s  flux  differencing.  In  the  simulations,  a  fifth-order 
WENO  scheme  is  used  on  the  vortex  tracking  meshes  to  study  the  advantages  of  using  a  higher- 
order  scheme  to  preserve  the  vortices  to  older  wake  ages  where  they  would  interact  with  the  ground. 
The  viscous  terms  are  computed  using  second-order  central  differencing.  The  Spalart-Allmaras 
turbulence  model  is  used  in  all  cases  for  RANS  closure.  Time  integration  is  performed  using  a 
second-order  implicit  backward-difference  method  scheme  with  dual  time-stepping.  The  Lower- 
Upper  Symmetric  Gauss-Siedel  (LUSGS)  method  is  used  for  the  inversion  for  the  sub-scale  rotor 
cases,  while  the  micro-scale  rotor  cases  use  the  diagonalized  approximate  factorization  framework. 


Results 

The  original  methodology  developed  in  the  current  task  has  shown  good  agreement  with  flow  mea¬ 
surements  made  for  a  small-scale  rotor  (referred  to  as  “micro-scale”)  hovering  close  to  the  ground. 
A  similar  application  of  RANS  methodology  using  S-A  turbulence  model  with  rotational  correc¬ 
tion  that  was  used  in  small-scale  rotor  results  in  excessive  smearing  of  the  flow  field  when  applied 
to  the  sub-scale  rotor.  The  excessive  smearing  for  the  sub-scale  rotor  flow  field  occurs  because  of 
high  turbulence  levels,  and  this  observation  can  be  made  by  comparing  the  eddy  viscosity  values 
for  the  small-scale  and  sub-scale  rotor  in  Figs.  1.15(a)  and  (b)  respectively.  The  S-A  turbulence 
model  seems  to  preserve  the  tip  vortices  well  for  the  small-scale  rotor  because  the  predicted  tur- 
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(a)  S-A  simulation  (b)  SA-DDES  simulation 

Figure  1.16:  Comparison  of  CFD  predicted  eddy  viscosity  contours  at  an  azimuthal  plane  of 
0°  for  the  sub-scale  rotor.  (Notice  the  different  contour  levels). 


bulence  levels  were  quite  low,  especially  in  the  regions  of  tip  vortex.  However,  turbulence  levels 
were  observed  to  be  an  order  of  magnitude  higher  in  the  sub-scale  rotor  flow  field. 

An  improvement  over  the  S-A  turbulence  modeling  is  made  by  utilizing  a  hybrid  RANS 
and  LES  methodology  and  one  such  approach  is  known  as  Delayed  Detached  Eddy  Simulation 
(DDES).  In  this  hybrid  approach  using  DDES,  one  can  model  the  small  scale  eddies  in  the  bound¬ 
ary  layer/near  wall  regions  and  can  resolve  the  larger  eddies  in  the  regions  of  massive  flow  separa¬ 
tion.  DDES  is  implemented  by  modifying  the  distance  function  in  S-A  model  that  is  chosen  from 
the  minimum  of  nearest  wall  distance  or  grid  spacing  for  each  grid  point.  The  grid  spacing  is  cho¬ 
sen  from  the  maximum  of  grid  spacing  in  all  three  directions  and  this  formulation  was  originally 
designed  for  isotropic  grids.  However,  the  grids  used  in  sub-scale  rotor  simulation  involve  large 
amount  of  difference  in  the  grid  spacing  in  the  three  directions  and  are  highly  anisotropic.  There¬ 
fore,  the  grid  spacing  is  based  on  a  function  that  is  dependent  on  the  aspect  ratio  of  grids.  This 
hybrid  methodology  is  referred  to  as  SA-DDES  in  the  following  results.  The  resulting  turbulence 
levels  with  SA-DDES  are  shown  in  Fig.  1.16(a)  and  are  compared  with  the  S-A  model  predicted 
levels  in  Fig.  1.16(b).  Notice  that  the  eddy  viscosity  levels  have  dropped  by  an  order  of  magnitude 
in  the  SA-DDES  simulations  compared  to  S-A  simulations. 

Figure  1.17  shows  a  comparison  of  the  CFD  (SA-DDES)  computed  vorticity  magnitude  values 
in  the  near  wake  region  for  different  tip  shapes  used  in  the  sub-scale  experiments.  It  can  be  seen 
that  the  tip  vortices  trailing  from  the  rectangular  and  swept  tips  look  similar  whereas  the  BERP- 
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(c)  BERP-like  (d)  Slotted 


Figure  1.17:  CFD  (SA-DDES)  predicted  tip  vortex  formation  from  different  blade  tip  shapes. 

like  tip  has  a  lower  vortex  strength.  The  slotted  tip  produces  a  highly  diffused  vortex  core,  which 
occurs  because  of  the  turbulence  produced  by  the  internal  slots  at  the  tips  of  the  blade  and  the 
subsequent  turbulent  mixing  produced  in  the  otherwise  laminar  tip  vortex. 

Figures  1 . 1 8(a)  and  (b)  compares  the  core  radius  growth  with  wake  age  for  the  four  tip  shapes  as 
measured  in  the  experiments  and  as  computed  by  CFD.  The  core  radius  was  non-dimensionalized 
by  the  rotor  blade  chord.  It  can  be  seen  that  the  core  radius  growth  agrees  well  for  all  blade  tip 
shapes.  The  initial  core  radius  for  the  BERP-like  tip  seems  to  be  higher  in  the  experiments  than 
the  computations  but  starts  to  compare  well  with  experiments  beyond  750°.  The  growth  of  core 
radius  for  the  slotted  tip  is  seen  to  be  experimentally  measured  only  up  until  a  wake  age  of  200°, 
however  CFD  can  still  track  the  core  radius  for  the  entire  range  of  wake  age  and  the  core  radius  is 
seen  to  double  at  800°  from  its  value  obtained  at  200°. 

Peak  swirl  velocity  were  compared  between  the  experimental  measurements  and  CFD  com- 
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(a)  Normalized  core  radius  comparison  (b)  Normalized  core  radius  comparison 

Figure  1.18:  Comparison  of  CFD  (SA-DDES)  predicted  core  radius  growth  as  a  function  of 
wake  age  for  different  tip  shapes. 


(a)  Non-dimensional  peak  swirl  velocity  (b)  Non-dimensional  peak  swirl  velocity 

Figure  1.19:  Comparison  of  CFD  (SA-DDES)  predicted  peak  swirl  velocity  as  a  function  of 
wake  age  for  different  tip  shapes. 
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putations  at  various  wake  ages,  as  shown  in  Figs.  1.19(a)  and  (b).  Both  the  experiments  and 
computations  showed  a  good  comparison  of  peak  swirl  velocity  until  360°  of  wake  age.  Beyond 
that,  the  effect  of  ground  leads  to  a  slight  increase  in  peak  swirl  velocity  in  experiments  whereas  it 
continues  to  drop  in  the  computations.  The  peak  swirl  velocity  variation  was  noticed  to  be  similar 
for  the  rectangular  and  swept  tip  shapes.  The  BERP-like  tip  shows  a  higher  peak  swirl  velocity 
at  early  wake  ages  in  computations  but  compares  well  at  older  wake  ages  with  the  experiments. 
The  slotted  tip  case  shows  higher  peak  swirl  velocity  in  computations  but  eventually  decays  and 
compares  well  with  the  experiments  after  400°  wake  age. 


(c)  BERP-like  (d)  Slotted 


Figure  1.20:  CFD  (SA-DDES)  predicted  pressure  coefficient  by  the  1-bladed  sub-scale  rotor 
blade  for  different  tip  shapes. 
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The  unsteady  pressure  variations  at  the  ground  plane  can  be  a  contributing  factor  for  sediment 
mobilization  during  brownout.  Figures  1.20(a)  and  (b)  show  pressure  coefficient  at  the  ground 
plane  for  the  all  tip  shapes.  The  flow  field  is  extracted  in  computations  from  an  instantaneous 
solution.  The  circular  region  extending  to  x/R  =  1 .0  of  red  contour  values  shows  the  high  pressure 
in  the  inboard  regions.  Beyond  that,  one  can  observe  lower  pressure  (blue  contour  levels)  because 
of  the  presence  of  tip  vortex  that  tends  to  induce  a  lower  pressure  region  as  the  tip  vortex  passes 
over  the  ground  plane.  Notice  the  striations/fluctuations  in  the  pressure  values  for  the  rectangular, 
swept  and  BERP-like  tip  shapes  which  are  not  present  in  the  slotted  tip. 
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Fuselage  Configuration  Effects  on  Rotorwash  and  Brownout 

Investigator(s):  Ganesh  Rajagopolan  (until  06/30/11) 
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Institution/Department:  University  of  Maryland,  Department  of  Aerospace  Engineering 

Graduate  Student(s):  Bharath  Govindarajan 

Contact  email(s):  leishman@umd.edu,  baeder@umd.edu 

Background  and  Technical  Challenges 

Brownout  is  degraded  visual  environment  (DVE)  problem  that  is  caused  by  clouds  of  dust  particles 
being  uplifted  by  a  rotorcraft  during  landings  and  takeoffs  in  desert-like  environments.  While 
the  motion  of  the  dust  particles  is  primarily  caused  by  the  effects  of  the  rotor  wake  vortices  and 
downwash  field  produced  by  the  rotor,  the  effects  of  the  airframe  shape  may  further  distort  the 
development  of  the  wake,  thereby  affecting  the  development  of  the  flow  at  the  ground  and  severity 
of  brownout  conditions.  Furthermore,  the  presence  of  the  tail  rotor,  which  also  produces  a  strong 
vortical  wake  system,  can  interact  with  both  the  airframe  flow  and  with  the  main  rotor  wake, 
thereby  potentially  affecting  the  development  of  brownout  conditions.  Although  mostly  based 
on  anecdotal  evidence  thus  far,  there  has  been  reason  to  suspect  that  the  shape  of  the  fuselage 
and/or  the  type  and  positioning  of  the  tail  rotor  has  at  least  some  measurable  influence  on  the 
development  brownout  conditions,  and  that  there  could  be  a  preferred  design  that  could  help  to 
minimize  brownout  severity,  e.g.,  to  allow  the  pilot  to  better  see  the  ground  through  the  dust  cloud. 

Technical  Approach 

The  objective  of  this  research  task  was  to  model  the  role  of  the  airframe  parts  besides  the  main  ro¬ 
tor  as  they  may  contribute  to  the  distortion  the  flow  around  a  rotorcraft  operating  in  ground  effect, 
and  to  better  understand  the  role  of  the  fuselage,  empennage,  tail  rotor,  etc.,  in  the  development 
brownout  and  DVE  conditions.  A  parallel  experimental  study  of  rotor/airframe/ground  effect  in¬ 
teractions  is  being  conducted  under  Tasks  1.1  and  2.2.  The  specific  objectives  of  this  task  are  to: 
1.  Develop  a  set  of  computationally  efficient  modeling  tools  to  study  the  aerodynamic  interaction 
of  the  main  rotor  system  when  combined  with  different  airframe/fuselage  configurations.  Con¬ 
sideration  is  being  given  to  modeling  the  flow  produced  by  tail  rotor,  as  well  as  multiple  rotor 
aircraft  such  as  tandems  and  coaxials;  2.  To  develop  modeling  tools  to  predict  and  study  parti- 
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cle  transport  around  complete  rotorcraft  configurations  that  are  operating  in  maneuvering  flight 
near  the  ground,  with  a  view  to  understanding  the  influence  of  the  airframe  and  tail  rotor  on  the 
development  of  brownout  conditions. 

An  initial  study  into  the  effects  of  fuselage  shape  on  the  groundwash  signature  below  a  rotor 
was  performed  by  Iowa  State  University  (ISU).  However,  the  ISU  approach  offered  limited  fidelity 
because  it  could  not  represent  the  vortical  features  of  the  rotor  wake  that  are  now  known  to  be 
critically  important  in  the  mobilization  and  uplift  of  dust  particles  (e.g.,  see  outcomes  from  Tasks 
1.1  and  Task  3.2).  In  the  ISU  approach,  the  inflow  to  the  rotor  disk  was  modeled  using  momentum 
source  terms  but  this  formulation  does  not  represent  the  vortices  produced  by  the  rotor  blades  and 
so  the  flow  near  the  ground  is  not  representative  of  the  real  physics  that  drives  the  brownout  prob¬ 
lem.  A  better  approach,  which  is  based  on  a  Lagrangian  free-vortex  rotor  wake  and  surface  singu¬ 
larity  (panel)  formulation  of  the  airframe,  has  now  been  developed.  This  latter  approach  is  more 
aligned  with  the  physics  of  predicting  wake  rotor  flows  in  ground  effect  and  the  onset  of  brownout 
conditions,  but  is  also  computationally  more  attractive  than  the  CFD  approach  (Task  1.2).  This 
approach  can  also  be  easily  interfaced  with  other  parallel  research  initiatives  of  this  MURI,  such 
as  the  ABATE  framework  described  under  Task  3.4,  as  well  as  the  various  other  mitigation  studies 
such  as  those  being  used  in  Task  3.3.  A  Lagrangian  free-vortex  methodology  was  also  employed 
to  model  the  tail  rotor  wake,  and  was  subsequently  coupled  with  the  main  rotor  and  surface  singu¬ 
larity  formulation  for  the  body.  This  foregoing  approach  has  provided  a  comprehensive  model  to 
study  the  particle  transport  physics  around  complete  rotorcraft  configurations,  and  has  been  used 
to  make  assessments  of  the  severity  of  the  brownout  cloud  as  well  as  to  investigate  possible  means 
of  brownout  mitigation. 

A  schematic  of  the  basic  methodology  is  shown  in  Fig.  1.21.  The  ground  plane  is  modeled  us¬ 
ing  an  image  system  to  ensure  that  flow  tangency  is  satisfied  at  all  points  over  the  ground  plane.  The 
Lagrangian  free-vortex  wake  methodology  (also  used  in  Task  3.2)  is  coupled  with  a  surface  singu¬ 
larity  methodology  made  up  of  a  series  of  quadrilateral  and/or  triangular  panels;  flow-tangency  is 
enforced  at  the  control  points  of  these  panels.  The  unsteady  pressure  on  the  surface  of  the  panel 
(i.e.,  at  the  control  point)  is  obtained  through  the  use  of  the  unsteady  Bernoulli’s  equation  (Kelvin’s 
equation),  where  the  unsteady  terms  represent  the  effects  associated  with  changes  in  panel  singu¬ 
larity  strength,  the  moving  vortices  in  the  near  wake,  and  the  rolled-up  tip  vortices  in  far  wake  of 
the  rotor.  In  some  case,  especially  when  vortices  come  near  the  panels,  the  unsteady  terms  can 
dominate  the  surface  loading.  The  present  approach  provides  a  computationally  attractive  predic¬ 
tion  of  the  vortical  rotor  wake,  the  wake  distortions  by  the  airframe,  and  the  resulting  unsteady 
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Main  rotor 


Figure  1.21:  Schematic  illustrating  the  modeling  approach  of  the  free-vortex  method  and 
coupled  rotor-fuselage  system  in  ground  effect  operation. 


three-dimensional  flow  as  it  develops  near  the  ground. 

The  modeling  of  the  tail  rotor  wake  has  also  been  performed  using  the  free-vortex  method.  The 
time-step  of  the  procedure,  At,  is  given  by  At  =  AT'/fl,  where  A*T  is  the  rotor  azimuthal  discretiza¬ 
tion  and  Q.  is  the  rotational  frequency  of  the  rotor.  Because  the  main  and  tail  rotor  rotate  at  different 
rotational  speeds,  the  coupling  of  information  (i.e.,  vortex  positions  and  strengths)  between  the  two 
rotors  can  be  potentially  realized  in  two  ways:  1.  The  azimuthal  discretization  of  the  two  rotors 
is  set  so  that  At  is  the  same  for  both  rotors,  which  is  called  tight  coupling  because  information  is 
passed  between  the  two  rotors  at  each  and  every  time  step;  2.  The  azimuthal  discretizations  are  set 
to  be  equal,  and  the  rotor  with  the  lower  At  marches  additional  time-steps  such  that  both  rotors  are 
now  at  the  same  physical  time  in  the  solution,  which  is  called  loose  coupling  because  information 
between  the  rotors  are  passed  at  predetermined  intervals.  A  schematic  of  this  latter  approach  is 
shown  in  Fig.  1.22.  The  loose  coupling  methodology  was  adopted  in  the  present  work  because 
it  was  found  to  be  computationally  faster  and  required  less  memory,  but  while  giving  equivalent 
accuracy  compared  to  the  tight  coupling  approach. 

The  rotor  airframe  model,  including  the  tail  rotor,  has  been  used  to  predict  flow  velocities  on 


42 


Q  Interpolated  solution  Diractson  of 

solulion  propagal'Dn 


Vertex  lilflmenl? 
Bitched  lo  blatters 
at  iJ>*pO 


On&  simutation  Limosiep 


■*- 


vonosr  (iiamonis 
attache  io  blades 

at 


Figure  1.22:  Schematic  illustrating  the  loose  coupling  methodology  used  to  compute  the  aero¬ 
dynamic  interactions  between  the  main  and  tail  rotors. 

the  ground  and  the  resulting  flow  solution  has  been  coupled  with  the  particle  uplift  framework  of 
Task  3.2  to  predict  the  formation  of  the  evolving  dust  clouds.  Parallelization  of  the  particle  uplift 
code  on  GPUs  (from  Task  2.6)  has  also  been  undertaken  so  as  to  reduce  the  computational  time 
while  also  preserving  the  accuracy  of  the  numerical  solution. 

The  trajectories  of  the  dust  particles  in  space  and  time  are  affected  not  only  by  the  aerodynamic 
influence  of  the  fuselage  and  vortex  filaments,  but  the  particles  may  impact  the  airframe  and  the 
blades,  and  the  resulting  collisions  may  significantly  alter  the  trajectories  of  the  particles.  Further¬ 
more,  the  presence  of  suspended  dust  particles  when  a  rotorcraft  operates  in  brownout  conditions 
not  only  limits  the  visibility  for  the  pilot  but  results  in  erosional  damage  of  the  rotor  blade,  which  is 
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Figure  1.23:  Validation  of  the  surface  singularity  method  with  the  time-averaged  pressure 
distribution  on  the  AGARD  fuselage  in  a  uniform  free-stream  flow. 

often  seen  on  the  metallic  leading  edge  and  the  composite  blade  skin  aft  of  the  leading  edge.  This 
type  of  damage  is  more  common  on  the  outboard  sections  of  the  blade,  where  the  local  flow  speeds 
relative  to  the  blade  are  highest.  This  erosive  damage  can  occur  rather  quickly,  and  may  render  the 
blade  aerodynamically  and  structurally  compromised  after  only  a  few  encounters  brownout  condi¬ 
tions.  The  particles  contained  in  the  flow  through  the  engine  can  also  erode  the  compressor  blades 
and  so  can  reduce  engine  performance.  Current  strategies  for  predicting  erosion  patterns,  and  the 
resulting  performance  degradation,  stems  primarily  from  wind-tunnel  (hot  erosion  and  cascade 
erosion)  testing,  coating  life  evaluations,  and  certain  types  of  particle  simulations.  Initial  work  has 
recently  been  conducted  in  this  task  to  develop  a  computational  framework  in  which  blade  abrasion 
effects  can  be  predicted  with  good  accuracy  and  at  relatively  low  computational  cost. 

Results 

Validation  of  the  Panel  Method 

Validation  of  the  surface  singularity  method  was  undertaken  using  simple  body  geometries 
such  as  a  sphere  (exact  solution)  and  also  the  slender  AGARD  body  (i.e.,  a  simplified  helicopter 
fuselage)  standard  test  case  for  which  there  are  measurements.  Figure  1.23  shows  the  variation 
of  pressure  coefficient  against  longitudinal  distance  along  the  AGARD  body  when  the  body  is 
operating  in  a  uniform  free-stream  flow  at  zero  angle  of  attack.  In  this  case,  the  rotor  blades 
were  removed  but  the  measurements  were  made  with  the  rotating  hub.  Each  color  represents  the 
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(a)  Retreating  side 


(b)  Top  surface 


Figure  1.24:  Validation  of  the  surface  singularity  method  for  the  AGARD  fuselage.  Steady 
pressure  distribution  on  the:  (a)  Retreating  side;  (b)  Top  surface. 


data  from  the  pressure  taps  installed  on  the  three  sides  of  the  body.  The  results  were  found  to  be 
in  reasonably  good  agreement  with  the  measured  data,  but  there  are  discrepancies.  In  this  case, 
regions  of  flow  separation  produced  by  the  rotor  hub  are  responsible  for  the  differences  shown. 

Figure  1.24  shows  the  time-averaged  pressures  on  the  top  and  retreating  side  of  the  AGARD 
fuselage  with  the  rotor  in  edgewise  forward  flight  at  an  advance  ratio,  /u,  of  0.05.  Again,  the  dif¬ 
ferences  in  the  predictions  and  measurements  on  the  top  surface  of  the  body  arise  because  of  the 
presence  of  the  rotor  hub  and  its  wake  was  not  modeled.  The  model  also  over-predicts  the  velocity 
at  the  front  end  of  the  fuselage  leading  to  decreased  pressure  values  on  the  retreating  side.  Over¬ 
all,  however,  the  results  are  in  reasonable  agreement  with  the  measurements  bearing  in  mind  the 
complexity  of  the  physical  problem  of  a  rotor  wake  interacting  with  a  solid  body. 


Interaction  of  Main  Rotor  and  Tail  Rotor 

Figures  1.25(a)  and  1.25(b)  show  the  longitudinal  and  lateral  views  of  the  predicted  wake 
geometry  with  and  without  the  presence  of  a  tail  rotor,  respectively,  for  the  hovering  rotor  operating 
out  of  ground  effect.  The  isolated  rotor  wake  was  noted  to  be  periodic  and  axisymmetric  for  nearly 
four  rotor  revolutions,  as  shown  in  Fig.  1.25(a).  At  later  wake  ages,  bundling  of  the  tip  vortices 
produced  aperiodic  wake  developments.  It  was  found  that  the  presence  of  the  tail  rotor  distorted 
the  main  rotor  wake  and  also  led  to  increased  levels  of  aperiodicity  in  the  wake  at  earlier  wake 
ages,  as  shown  in  Fig.  1.25(b). 

The  tail  rotor  vortex  filaments  (highlighted  in  blue)  were  entrained  through  the  downwash  of 
the  main  rotor  where  it  interacted  with  the  tip  vortex  filaments  of  the  main  rotor  at  very  early  wake 
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(b)  SMR  with  TR 


Figure  1.25:  Longitudinal  and  lateral  wake  geometries  for  (a)  an  isolated  SMR  and  (b)  SMR 
with  a  TR  operating  in  hover  out  of  ground  effect. 

ages.  As  compared  to  the  tip  vortex  structure  of  the  isolated  rotor,  the  presence  of  the  tail  rotor  in¬ 
troduced  mild  disturbances  in  the  early  ages  of  the  main  rotor  wake  developments,  which  resulted 
in  the  bundling  of  main  rotor  tip  vortices  as  early  as  the  first  rotor  revolution.  The  induced  velocity 
of  the  tail  rotor  on  the  main  rotor  wake  also  resulted  in  a  lateral  wake  skew  tilt  exhibited  by  the 
wake  of  the  main  rotor,  as  shown  in  Fig.  1.25(b). 


Interaction  of  Main  Rotor  and  Fuselage 

Two  representative  fuselage  shapes  were  considered  to  approximately  quantify  the  effect  of  the 
fuselage  shape  on  the  development  of  the  main  rotor  wake  and  the  effects  produced  at  the  ground: 
1.  A  geometrically  simplified  helicopter  fuselage,  referred  to  as  the  UMD-AGARD  fuselage,  and 
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Figure  1.26:  Panel  representations  of  the  two  fuselage  shapes. 


(a)  Isolated  SMR 


(b)  UMD-AGARD  fuselage 


Tip- vortex  filaments 
perturbed  by  the 
presence  of  fuselage 


(c)  Cylindrical  fuselage 


Figure  1.27:  Comparison  of  the  rotor  wake  geometries  in  ground  effect  operation  for  the:  (a) 
isolated  SMR,  (b)  SMR  with  the  UMD-AGARD  body,  and  (c)  SMR  with  the  cylindrical  body. 

2.  A  fuselage  with  a  cylindrical  cross-section  and  hemispherical  caps,  as  shown  in  Fig  1.26.  While 
the  cylindrical  body  was  of  uniform  cross-section,  the  UMD-AGARD  fuselage  had  a  cylindrical 
main  body  that  tapered  at  the  nose  and  a  cylindrical  tail  boom  body  that  tapered  at  the  tail;  see 
Fig.  1.26. 

Figure  1.27  shows  the  wake  geometries  of  the  isolated  SMR,  as  well  as  the  SMR  in  the  pres¬ 
ence  of  the  two  fuselage  shapes.  The  results  are  shown  in  terms  of  the  filament  locations  as  they 
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Figure  1.28:  Variation  of  the  figure  of  merit  of  the  rotor  with  and  without  the  presence  of  a 
fuselage. 

intersect  the  longitudinal  plane.  For  the  UMD-AGARD  fuselage,  the  trailed  tip  vortices  were  only 
minimally  distorted  as  they  slowed  and  stretched  around  the  tail,  but  this  process  clearly  initiates  a 
wake  disturbance  and  it  can  be  seen  that  these  disturbances  are  also  propagated  to  the  front  of  the 
rotor  wake.  This  means  that  local  wake/fuselage  interactions  may  have  some  effects  on  the  overall 
flow  environment  at  the  ground. 

A  study  was  also  conducted  to  quantify  the  influence  of  the  fuselage  on  the  performance  of 
the  rotor.  Figure  1.28  shows  the  variation  of  the  rotor  figure  of  merit  against  the  blade  loading 
coefficient,  Ct/g,  with  and  without  the  AGARD  fuselage  (see  Fig.  1.27(b))  operating  in  hover  out 
of  ground  effect.  The  solid  circles  indicate  experimental  measurements,  while  the  dashed  lines 
represent  FVM  predictions.  As  the  rotor  thrust  increases,  the  presence  of  the  fuselage  produces 
an  increase  in  the  rotor  figure  of  merit;  an  increase  of  about  5%  at  Cj/g  =  0. 1  is  clearly  significant. 

Brownout  Dust  Cloud  Simulations  -  Effects  of  Fuselage 

The  nominal  landing  profile  for  the  simulations  used  is  shown  in  Fig.  1.29.  The  maneuver 
begins  with  the  rotorcraft  in  steady  forward  flight  at  50  kts  and  at  a  descent  rate  of  100  ft  min-1. 
The  helicopter  then  executes  a  flare  maneuver  and  decelerates  into  a  steady  hover  with  the  rotor 
plane  stabilizing  at  a  height  of  one  rotor  radius  above  the  ground.  The  duration  of  the  entire 
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Figure  1.29:  Approach  profile  for  landing:  (a)  Rotor  hub  center  positions;  (b)  Shaft  tilt 
angles;  (c)  Overall  maneuver. 

approach  is  about  35  s,  and  time  is  referenced  to  the  start  of  the  brownout  simulation,  as  shown  in 
Fig.  1.29. 

Figure  1.30  shows  predictions  of  the  dust  clouds  as  produced  by  the  isolated  SMR  at  the  end 
of  the  simulation,  for  the  SMR  with  the  UMD-AGARD  fuselage,  and  for  the  SMR  with  the  cylin¬ 
drical  fuselage,  respectively.  As  shown  in  Fig.  1.30,  the  dust  cloud  can  be  decomposed  into  two 
distinct  regions,  an  outer  structure  and  an  inner  structure.  The  outer  structure,  highlighted  in  red, 
develops  during  the  formation  and  propagation  of  the  ground  vortex.  These  structures  contain  par¬ 
ticles  that  were  uplifted  through  the  induced  velocity  field  produced  by  the  propagating  vortex  and 
also  by  those  particles  that  were  trapped  within  the  ground  vortex  and  so  ejecting  more  particles  by 
bombardment  mechanisms.  The  outer  structure  of  the  cloud  appears  to  be  sensitive  to  what  hap¬ 
pens  during  the  transitional  phase  from  forward  flight  to  hover,  in  which  wake  distortions  caused 
by  the  presence  of  the  fuselage  translate  into  more  global  flow  distortions  and  further  intricacies  in 
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(b)  UMD-AGARD  fuselage 


(c)  Cylindrical  fuselage 


Figure  1.30:  Simulated  dust  cloud  profiles  for  the  (a)  Isolated  rotor;  (b)  Rotor  and  UMD- 
AGARD  fuselage;  and  (c)  Rotor  and  cylindrical  fuselage. 

the  structure  of  the  dust  cloud. 

The  inner  structure  of  the  cloud,  highlighted  in  blue  in  Fig.  1.30,  is  characteristic  of  the  dome¬ 
shaped  dust  cloud  that  based  on  field  tests  in  brownout  conditions  is  often  seen  to  engulf  the  rotor. 
Because  the  presence  of  fuselage  only  slightly  modifies  the  rotor  wake  in  hover  (see  Fig.  1.27),  the 
inner  structure  of  the  cloud  is  seen  to  be  fairly  similar  in  all  cases.  However,  it  is  recognized  that 
these  solutions  may  change  when  flow  separation  effects  produced  by  the  fuselage  are  included 
into  the  modeling,  which  cannot  be  predicted  in  this  case  because  of  the  limitations  of  the  potential 
flow  model. 


Brownout  Dust  Cloud  Simulations  -  Effects  of  Rotor  Configurations 

Figure  1.31  shows  the  dust  cloud  structures  formed  by  the  SMR  (with  TR),  coaxial,  tandem, 
and  side-by-side  rotor  configurations  at  the  end  of  the  simulated  landing  maneuver.  The  dust  cloud 
produced  by  the  SMR  shown  in  Fig.  1.31(a)  is  similar  in  structure  to  that  of  the  isolated  SMR,  but 
with  minor  differences.  The  aerodynamic  influence  of  the  tail  rotor  (shown  in  Fig.  1.25)  causes 
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Figure  1.31:  Dust  clouds  formed  by  a  SMR  (with  TR),  coaxial,  tandem,  and  side-by-side 
rotors  at  the  end  of  the  simulated  landing  maneuvers. 

bundling  of  the  tip  vortices  and  makes  the  resulting  dust  cloud  considerably  less  axisymmetric. 
The  inclusion  of  the  tail  rotor  also  modifies  the  dust  cloud  by  causing  airborne  particles  to  be 
entrained  into  its  own  downwash  field. 

Figure  1.31(b)  shows  the  dust  cloud  produced  by  the  coaxial  rotor  system.  Unlike  other  rotor 
configurations,  the  lower  rotor  of  the  coaxial  significantly  increases  the  overall  slipstream  veloci¬ 
ties  in  the  wake,  which  causes  the  tip  vortices  to  impinge  on  the  ground  earlier  during  the  maneuver 
than  was  seen  for  the  other  rotor  configurations.  The  vortex  filaments  from  the  upper  and  lower 
rotor  bundle  and  convect  along  the  ground,  uplifting  large  quantities  of  dust  through  bombardment 
mechanisms.  The  coaxial  rotor  configuration  produces  a  dome  shaped  dust  cloud  that  is  more 
axisymmetric  compared  to  the  other  configurations. 

The  dust  cloud  produced  by  the  side-by-side  rotor  configuration  is  shown  in  Fig.  1.31(d).  It 
is  evident  that  this  configuration  uplifts  large  quantities  of  dust  that  is  distributed  on  all  sides,  po- 
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tentially  causing  severe  visual  obscuration  all  around  the  aircraft.  As  was  seen  with  the  tandem 
rotor  configuration,  a  strong  asymmetry  in  the  development  of  the  dust  cloud  exists  because  of  the 
positioning  of  the  rotors.  However,  in  this  case,  the  tip  vortex  filaments  from  each  rotor  bundle 
more  readily  and  more  frequently  and  so  forming  a  much  stronger  ground  vortex,  which  in  this 
case  convects  towards  the  rear  of  the  rotorcraft,  uplifting  large  quantities  of  dust  in  the  process. 
The  side-by-side  configuration  was  the  only  configuration  where  large  quantities  of  dust  were  con¬ 
tinuously  mobilized  in  hover,  which  is  because  of  the  relatively  higher  flow  velocities  near  formed 
the  ground  and  the  continuous  bundling  of  the  tip  vortex  filaments  in  the  outer  regions  of  the  wake. 

Particle  Collisions  with  the  Blade 

Blade  erosion  occurs  because  of  the  interaction  of  the  rotor  blades  with  the  suspended  partic¬ 
ulates  causing  wear  and  tear  of  the  leading  edge  of  the  blades.  In  a  brownout  environment,  this 
erosion  process  is  greatly  accelerated.  To  prevent  damage  from  dust/sand  erosion,  metallic  abra¬ 
sion  strips,  typically  composed  of  stainless  steel  or  titanium,  are  bonded  to  the  blade  leading  edge 
to  serve  as  a  hard  surface  that  absorbs  the  kinetic  energy  of  the  sand  particle.  These  impacting 
sand  particles,  however,  slowly  erode  away  the  metallic  material.  This  problem  is  a  bigger  issue 
near  the  outboard  blade  tip,  where  the  high  blade  velocity  means  that  the  dust  particles  impact  with 
greater  kinetic  energy.  Replacement  of  these  damaged  abrasion  strips  and  caps  on  the  rotor  blades, 
if  possible,  leads  to  increased  maintenance  cost  and  reduced  operational  readiness. 

For  example,  Figs.  1.32(a)-l. 32(d)  shows  the  trajectory  of  particles  of  sizes  10  /./m,  20  /./m, 
50  /./m,  and  70  /um,  respectively,  as  they  convect  around  a  NACA  0015  airfoil.  The  airfoil  was 
operating  at  10°  angle  of  attack.  The  time  taken  by  a  spherical  particle  to  respond  to  changes  in  the 
flow,  i.e.,  particle  time  constant,  is  proportional  to  the  square  of  particle  diameter  djr  Therefore, 
larger  the  particle  time  constant,  longer  it  takes  to  alter  its  trajectory,  and  consequently  it  has 
a  higher  probability  of  colliding  with  the  airfoil.  Figure  1.32(b)  shows  that  under  certain  flow 
conditions,  multiple  collisions  of  the  same  particle  with  the  airfoil  surface  can  occur  (trajectory 
highlighted  in  black). 
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Figure  1.32:  Trajectories  of  particles  colliding  with  a  NACA  0015  airfoil  at  10°  for  particles 
sizes  of  (a)  10  ji/m,  (b)  20  ji/m,  (c)  50  /urn,  and  (d)  70  ^um. 
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Part  2:  Fundamentals  of  Particle  Suspension 

•  Task  2.1  -  Non-Uniform,  Near-Bed  Flow  Field  with  Impinging  Rotorwash 

•  Task  2.2  -  Two-Phase  PIV  Near- Wall  Bounded  Turbulent  Flows 

•  Task  2.3  -  Fundamental  Two-Phase  Measurements  in  Brownout  Fluid  Mechanics 

•  Task  2.4  -  Large-Eddy  Simulation  of  the  Interaction  Between  Vortices  and  the  Ground 

•  Task  2.5  -  Two-Phase  Large  Eddy  Simulation  Using  the  Mesoscopic  Eulerian  Formalism 

•  Task  2.6:  Effects  of  Model  Scaling  on  Sediment  Transport 

•  Task  2.7:  Turbulent  Near-Bed  Flows  Under  Impinging  Vorticity 
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Task  2.1 


Non-Uniform,  Near-Bed  Flow  Field  Associated  With  Impinging  Rotorwash 

Investigator(s):  W.  Brian  Dade,  Benoit  Cushman  Roisin 
Institution/Department:  Dartmouth  College/Earth  Sciences 
Contact  email(s):  W.B.Dade@dartmouth.edu 

The  work  initially  being  performed  under  this  task  is  now  being  continued  under  two  other 
tasks,  namely  Task  2.6  (Effects  of  Model  Scaling  on  Sediment  Transport)  and  Task  2.7  (Turbulent 
Near-Bed  Flows  Under  Impinging  Vorticity).  Refer  to  these  tasks  for  further  information. 
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Task  2.2 


Dual-Phase  PIV  in  Two-Phase  Near- Wall  Bounded  Turbulent  Flows 

Investigator(s):  J.  Gordon  Leishman 

Institution/Department:  University  of  Maryland/Department  of  Aerospace  Engineering 
Graduate  Student(s):  Anish  Sydney,  Jaime  Reel,  Nathan  Doane 
Contact  email(s):  leishman@umd.edu 

Background  and  Technical  Challenges 

A  helicopter  rotor  wake  comprises  a  highly  three-dimensional  flow  with  locally  high  unsteady 
effects  produced  by  the  tip  vortices  from  the  rotating  blades.  Understanding  this  complicated 
vortically-dominated  rotor  wake  structure  becomes  even  more  difficult  as  it  approaches  a  ground 
plane.  Complexities  arise  because  the  vortices  undergo  stretching  and  bundling  as  they  interact 
with  the  developing  wall  jet  like  flow  along  the  ground.  The  various  effects  produced  by  the  tip 
vortices,  which  include  locally  high  shear  stresses  and  pressure  forces,  are  then  responsible  for 
mobilizing  and  uplifting  loose  sediment  particles.  The  mechanisms  by  which  sediment  particles 
are  uplifted  are  summarized  in  Fig.  2.1. 

The  ability  to  properly  understand  the  resulting  two-phase  flow  field,  as  well  as  how  it  is  af¬ 
fected  by  various  geometric  and  operating  parameters  of  the  rotor,  is  fundamental  to  understanding 
the  brownout  problem.  To  this  end,  the  objectives  of  this  task  are  to  develop  and  implement  di¬ 
agnostic  techniques  to  measure  the  fundamental  two-phase  flow  physics  underlying  the  brownout 
phenomenon.  The  main  goals  are  to  quantify  the  interactions  of  the  rotor  wake  with  a  ground  plane 
covered  with  a  mobile  sediment  bed,  determine  what  particular  characteristics  of  the  resulting  flow 
are  responsible  for  sediment  mobility  and  uplift,  and  then  how  the  resulting  dual-phase  flow  field 
is  altered  by  aspects  of  the  rotor  geometry  and  its  operating  condition. 

Technical  Approach 

Flow  field  measurements  have  been  conducted  in  the  controlled  environment  of  a  dust  chamber 
using  a  prototypical  flow  generated  by  a  small  rotor  system.  A  schematic  of  the  experimental 
setup  being  used  in  this  task  is  shown  in  Fig.  2.2.  The  rotor  has  a  radius  of  85  mm  (3.346  in), 
the  blades  having  a  cambered  plate  section.  Particle  image  velocimetry  (PIV)  experiments  were 
conducted  using  both  one-bladed  and  two-bladed  rotor  configurations.  The  teetering  hub  design 
allowed  the  second  blade  to  be  interchanged  with  an  equivalent  counter-mass  to  obtain  a  one- 
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Figure  2.1:  Schematic  showing  the  different  modes  of  sediment  mobilization  and  uplift  below 
a  rotor. 


bladed  rotor.  An  adjustable  height  ground  plane  was  arranged  so  that  it  was  parallel  to  the  rotor 
tip-path-plane.  Several  ground  planes  have  been  used  in  this  work,  including  one  fitted  with  a  row 
of  pressure  transducers  to  allow  measurements  of  the  unsteady  pressures  over  the  ground.  These 
transducers  were  spaced  1.3  cm  apart  except  between  r/R  =  1-2  where  they  were  spaced  at  0.7  cm 
to  obtain  a  higher  measurement  resolution.  Experiments  were  performed  with  the  rotor  hovering 
at  one  rotor  radius  above  the  ground  plane. 

Measurements  were  initially  made  with  the  laser  sheet  oriented  in  a  vertical  plane,  as  shown  in 
Fig.  2.2.  However,  to  better  understand  the  three-dimensionality  of  the  flow,  the  laser  light  sheet 
was  aligned  to  be  parallel  to  the  ground  plane  and  the  flow  imaged  by  viewing  downward  onto  the 
laser  sheet.  These  regions  were  also  used  to  study  the  differences  in  the  wall-parallel  velocities  near 
the  ground  plane  that  were  caused  by  the  presence  of  body  shapes  (representing  airframes)  in  the 
rotor  wake.  The  bodies  had  cross-sectional  shapes  that  were  circular,  elliptical,  and  rectangular, 
respectively,  as  shown  in  Fig.  2.3.  These  bodies  are  smaller  versions  of  those  used  in  Task  1.1. 
Dual-phase  experiments  were  also  performed  in  both  the  vertical  and  horizontal  planes. 

To  date,  most  studies  have  consisted  of  experiments  with  the  rotor  operating  at  various  fixed 
heights  above  the  ground.  In  practice,  however,  rotorcraft  encounter  brownout  conditions  during 
takeoffs  and  landings,  which  are  dynamic  flight  maneuvers.  While  certain  experiments  have  been 
performed  at  other  institutions  to  examine  the  problem  of  particle  pickup  induced  by  a  rotor  oper¬ 
ating  in  a  wind  tunnel,  the  wind  tunnel  approach  cannot  properly  simulate  the  boundary  conditions 
at  the  ground  plane.  Therefore,  a  rig  was  developed  in  the  present  work  that  allows  the  rotor  to 
be  translated  relative  to  the  ground  plane;  see  Fig.  2.4.  This  rig  allows  for  further  validation  of 
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Figure  2.2:  Schematic  showing  the  experimental  setup  used  to  perform  the  dual-phase  flow 
measurements. 


Figure  2.3:  Cross-sectional  views  of  the  body  shapes:  (a)  circular,  (b)  elliptical,  and  (c)  rect¬ 
angular. 

the  modeling  outcomes  (Tasks  3.2  and  3.4),  providing  the  necessary  range  of  motion  to  perform  a 
variety  of  simulated  flight  profiles. 

A  range  of  rotor  operating  conditions  and  heights  above  the  ground  have  been  used  to  expose 
the  details  of  the  wake  interactions  with  the  ground,  and  also  to  isolate  the  different  sediment  uplift 
mechanisms  that  ultimately  may  affect  the  formation  of  the  dust  cloud.  Dual-phase  particle  image 
velocimetry  (PIV)  was  used  to  measure  the  resulting  flows.  The  mobilized  sediment  particles  were 
identified  and  tracked  through  the  flow  by  using  particle  tracking  velocimetry  (PTV). 

Results 

The  mechanisms  previously  identified  in  Fig.  2.1  cause  loose  particles  on  the  bed  to  be  mobilized 
and  transported  into  the  rotor  flow  field.  The  particles  in  the  near-wall  region  experience  several 
types  of  forces,  including  shear,  pressure,  inter-particle,  and  gravitational.  The  pronounced  vertical 
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Figure  2.4:  Image  of  the  two-axis  traverse  system  used  to  simulate  descent  and  forward  flight. 

uplift  of  particles  that  was  observed  in  a  region  strongly  affected  by  the  predominantly  wall-parallel 
velocities  and  associated  boundary  layer  shear  suggests  additional  body  forces  could  be  acting  on 
the  particles.  For  example,  the  forces  on  the  particles  resulting  from  the  low  pressures  and  pressure 
gradients  near  the  vortex  core  may  be  responsible  for  the  more  wall-normal  particle  trajectories  that 
have  been  observed. 

To  quantify  the  pressure  induced  on  the  ground  plane,  experiments  were  performed  with  the 
rotor  operating  above  the  ground  plane  fitted  with  pressure  transducers.  A  representative  pressure 
signal  from  a  transducer  located  at  r/R  =  1.2  is  shown  in  Fig.  2.6.  This  radial  location  is  near  where 
the  vortices  impinge  on  the  ground  and,  therefore,  shows  a  more  periodic  signal,  i.e.,  before  the 
vortices  have  begun  to  diffuse  and/or  interact  with  each  other.  Notice  that  there  is  a  sharp  peak  in 
pressure  approximately  every  blade  passage,  which  is  caused  by  the  close  passage  of  a  tip  vortex 
over  a  pressure  transducer. 

To  further  examine  the  pressure  variation  along  the  ground,  Fig.  2.7  shows  pressure  signals 
from  all  of  the  sensors  as  a  function  of  wake  age.  Notice  that  the  largest  fluctuating  pressures  are 
seen  closer  to  the  rotor  ( r/R  =  0.5),  where  there  is  more  vertical  downwash  flow  from  the  rotor. 
The  pressure  decreases  as  the  blade  tip  vortices  approach  the  ground  plane  and  begin  to  stretch  and 
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Figure  2.5:  The  dust  chamber  in  which  the  flow  experiments  are  performed. 

expand  radially  outward.  In  the  region  of  vortex  impingement  on  the  ground,  the  pressure  become 
strongly  negative  (i.e.,  a  suction  pressure)  before  settling  back  to  near  ambient  conditions  further 
downstream.  However,  fluctuations  were  still  seen  for  all  radial  locations.  Near  the  rotor  these 
fluctuations  were  strongest,  but  the  effects  of  vortical  flows  and  turbulence  were  noted  to  cause 
pressure  variations  as  far  downstream  as  r/R  =  4. 

The  three-dimensional  nature  of  the  flow  field  has  been  examined  in  a  horizontal  plane  at  z/R  = 
0.1  (just  above  the  ground  plane),  the  results  being  shown  in  Fig.  2.8  as  contours  of  total  velocity. 
The  tip  vortices  manifest  as  localized  regions  of  higher  flow  velocities,  such  as  the  red  region 
marked  as  Filament  2  in  Fig.  2.8.  The  significant  three-dimensionality  of  the  rotor  wake  flow  at 
the  ground  is  evident.  At  x/R  =  0.6-0.8,  a  lower  velocity  region  under  the  rotor  was  observed,  this 
region  being  just  upstream  of  the  point  where  the  rotor  wake  boundary  had  begun  to  impinge  on 
the  ground.  By  x/R  =  0.9,  however,  the  rotor  wake  had  turned  radially  outward.  Just  outboard 
of  this  radial  location,  two  vortex  filaments  were  observed  to  impinge  on  the  ground,  which  are 
marked  as  Filament  1  and  Filament  2,  respectively.  Notice  that  the  impingement  of  these  vortices  is 
fairly  radially  asymmetric  along  the  length  of  the  filaments.  Beyond  x/R  =  1.5,  the  vortex  filaments 
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Figure  2.6:  Variation  of  surface  pressure  produced  at  r//?  =  1.2  by  a  one-bladed  rotor  hover¬ 
ing  at  z/R  =  1.0. 


Figure  2.7:  Time-histories  of  surface  pressures  produced  by  a  one-bladed  rotor  hovering  at 
z/R  =  1.0. 

were  not  distinguishable  from  each  other  because  they  had  either  merged  together  or  had  diffused, 
causing  a  highly  unsteady,  radially  asymmetric  flow  environment. 

The  effects  of  placing  a  body  (or  airframe)  in  the  rotor  wake  on  the  mobilization  and  uplift 
of  sediment  has  been  examined.  Figure  2.9  shows  the  dual-phase  flow  produced  by  the  isolated 
rotor  compared  to  the  flow  with  the  rectangular  body  present.  The  background  contours  are  of 
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Figure  2.8:  Total  velocity  contours  for  the  isolated  rotor  hovering  at  z/R  =  1.0  in  a  horizontal 
plane  at  z/R  =  0.1. 


total  velocity  and  the  sediment  particle  locations  are  identified  as  white  circles;  the  circles  have 
been  exaggerated  in  size  for  clarity.  Notice  that  these  results  are  in  a  vertical  plane  coincident  with 
the  longitudinal  centerline  of  the  bodies,  and  show  the  tail  region  of  the  rectangular  body.  The 
carrier-phase  measurements  (i.e.,  the  rotor  wake  itself)  shows  that  the  body  significantly  altered 
the  flow  at  the  tail  section.  The  tip  vortices  interacted  with  the  body  at  relatively  young  wake  ages, 
directly  impinging  on  the  body  surface  and  causing  the  vortices  to  rapidly  diffuse  in  some  cases. 
However,  relatively  high  regions  of  velocity  near  the  ground  suggest  that  some  remaining  vorticity 
does  persist  to  older  wake  ages.  In  general,  the  vortical  flow  at  the  ground  with  a  body  present 
was  found  to  be  significantly  more  diffused  compared  to  that  for  the  isolated  rotor.  Of  particular 
interest,  however,  is  that  the  flow  field  produced  by  the  presence  of  the  body  appears  to  uplift  more 
particles  than  for  the  isolated  rotor  flow. 

While  somewhat  counterintuitive,  this  latter  observation  is  a  consequence  of  the  higher  peak 
velocity  that  was  produced  at  the  ground  by  the  isolated  rotor.  These  higher  flow  velocities  con¬ 
verted  the  particles  more  radially  outward  away  from  the  rotor  before  they  could  be  picked  up  and 
suspended  by  the  blade  tip  vortices.  With  the  bodies  present,  however,  the  near-wall  flow  reached 
a  lower  peak  velocity  at  all  radial  locations  away  from  the  rotor.  While  still  above  the  threshold 
conditions  required  to  mobilize  the  particles,  the  near-wall  flow  did  not  convert  the  particles  radi¬ 
ally  outward  as  rapidly  or  as  far  as  with  the  isolated  rotor.  Instead,  the  vortical  flow  near  the  wall 
(which  actually  contained  less  coherent  vorticity  than  was  produced  by  the  isolated  rotor)  resulted 
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(a)  Dual-phase  measurements  for  the  isolated  rotor 
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(b)  Dual-phase  measurements  for  the  rectangular  body 

Figure  2.9:  Instantaneous  dual-phase  measurements  for  the  isolated  rotor  with  the  rectan¬ 
gular  body. 

in  more  particles  being  uplifted  closer  to  the  rotor. 

The  measurements  made  in  the  vertical  imaging  plane  showed  that  the  flow  in  the  tail  region  of 
the  body  was  significantly  distorted.  Therefore,  the  flow  field  was  measured  over  several  horizontal 
planes  to  better  quantify  such  flow  distortions;  the  single-phase  results  for  z/R  =  0.1  are  shown  in 
Fig.  2.10  for  the  isolated  rotor  and  for  the  rotor  with  the  rectangular  body,  respectively.  The 
background  contours  are  of  total  velocity  time-averaged  over  1,000  contiguous  PIV  realizations 
of  the  flow.  The  three-dimensional  flow  in  this  region  was  seen  to  be  significantly  altered  by  the 
presence  of  the  bodies,  causing  lower  average  flow  velocities  and  more  radial  flow  asymmetries 
compared  to  those  produced  by  the  isolated  rotor.  Another  interesting  difference  from  the  isolated 
rotor  flow  is  the  formation  of  a  localized  flow  adjacent  to  the  tail  of  the  body,  which  is  identified  in 
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(a)  Isolated  rotor  (b)  Rectangular  body 


Figure  2.10:  Time-averaged  total  velocity  contours  for  the  isolated  rotor  and  below  the  tail 
of  each  body  at  z/R  =  0.1. 


x/R 

(a)  Particle  concentration  -  isolated  rotor 


x/R 

(b)  Particle  concentration  -  rectangular  body 


Figure  2.11:  Time-averaged  particle  concentration  at  z/R  =  0.3  for  the  isolated  rotor  and  at 
the  tail  of  the  rectangular  body. 


Fig.  2.10(b). 

To  examine  the  three-dimensional  nature  of  the  the  particle  field  with  the  bodies  present,  dual¬ 
phase  tests  were  conducted  in  a  horizontal  plane  at  z/R  =  0.3.  The  particle  locations  were  binned 
and  time-averaged  over  1,000  contiguous  realizations  of  the  flow  to  produce  concentration  maps; 
see  Fig.  2.11.  The  particle  concentrations  produced  by  the  isolated  rotor  indicated  that  particles 
were  being  uplifted  and  suspended  at  x/R  =  2.4  and  then  continuing  downstream.  With  the  body 
present,  however,  the  particles  were  suspended  closer  to  the  tail  region  at  x/R  ~1.0.  The  closer  sus¬ 
pension  of  particles  occured,  in  part,  because  the  localized  flow  near  the  mobilized  large  quantities 
of  particles  relatively  close  to  the  rotor.  In  general,  it  was  found  that  the  flows  produced  with  the 
bodies  present  caused  non-uniform  particle  concentrations  that  were  formed  closer  to  the  rotor. 
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Figure  2.12:  Two-bladed  rotor  in  level  forward  flight  at  /u  «  0.03. 


In  practice,  rotorcraft  are  generally  performing  take-off  or  landing  maneuvers  when  they  en¬ 
counter  brownout  conditions.  Several  experiments  have  recently  been  conducted  at  a  variety  of 
simulated  forward  flight  conditions  using  a  rotor  traverse  rig.  It  was  found  that  in  low-speed  for¬ 
ward  flight  in  ground  effect  the  flow  produced  by  the  rotor  goes  through  three  stages:  1.  The  wake 
was  first  convected  out  in  front  of  the  rotor;  2.  The  rotor  then  approaches  a  translational  speed  at 
which  the  flow  begins  to  recirculate  back  into  the  rotor  disk  causing  a  ground  vortex  to  form;  3. 
The  rotor  then  attains  enough  forward  speed  to  sweep  the  wake  back  behind  the  rotor. 

A  very  low-speed  case  at  n  «  0.03  is  shown  in  Fig.  2.12.  The  results  show  that  some  amount 
of  recirculation  was  present  in  the  flow,  but  the  wake  remained  close  to  the  the  ground  to  form  a 
wall  jet,  as  was  seen  in  hover.  A  case  for  jj  ~  0.07  is  shown  in  Fig.  2.13.  At  this  slightly  higher 
forward  speed  the  results  show  that  the  rotor  wake  detached  from  the  ground  and  was  recirculated 
back  into  the  rotor  disk.  Finally,  for  /u  ~  0.09,  the  rotor  wake  and  its  vortices  were  swept  well  back 
behind  the  rotor;  see  Fig.  2.14. 

A  ground  vortex  was  seen  to  develop  in  front  of  the  rotor  before  the  wake  was  swept  back  so 
a  smaller  region  of  interest  was  used  to  obtain  higher  resolution  flow  measurements.  As  shown 
in  Fig.  2.15(a),  contours  of  velocity  magnitude  suggest  that  rather  than  forming  a  typical  wall-jet 
like  flow  seen  in  hover,  the  flow  was  reingested  through  the  rotor  disk.  In  the  vorticity  plot  shown 
in  Fig.  2.15(b),  the  recirculating  flow  has  developed  positive  vorticity.  Furthermore,  the  vortical 
structures  persisted  to  old  enough  wake  ages  that  they  could  be  identified  even  as  the  flow  had  been 
turned  upward  towards  the  rotor  disk. 
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Figure  2.13:  Two-bladed  rotor  in  level  forward  flight  at  p  «  0.07. 
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Figure  2.14:  Two-bladed  rotor  in  level  forward  flight  at  p  ta  0.09. 
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Figure  2.15:  Velocity  magnitude  and  vorticity  of  ground  vortex  ahead  of  the  rotor. 
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Background  and  Technical  Challenges 

To  make  progress  in  the  prediction  of  particle  suspension  and  sedimentation  within  coupled  particle¬ 
laden  flows  relevant  to  the  problem  of  rotorcraft  brownout,  detailed  characterization  of  the  micro¬ 
scale  mechanics  is  needed  within  a  prototypical  flow  that  captures  the  essence  of  the  rotorcraft 
wake/ground  interactions.  Detailed  dynamic  characterization  of  the  local  particle  concentration 
and  velocity  statistics  are  required  relative  to  the  large-scale  turbulent  vortical  structures  that  are 
inherent  characteristics  of  brownout  flow  fields.  Such  measures  can  be  used  to  obtain  local  particle 
flux  in  response  to  turbulent  fluid  stresses,  while  further  details  on  the  coupling  can  be  derived  from 
measures  of  the  particle  to  fluid  velocity  cross-correlations.  The  challenges  involved  in  obtaining 
this  type  of  information  are  numerous:  1.  Well-controlled  particle  characteristics  are  needed  to 
generate  repeatability  in  the  experiments,  as  the  suspension  characteristics  are  known  to  be  sensi¬ 
tive  to  shape,  surface  roughness,  etc.;  2.  Making  reliable  high-resolution  simultaneous  two-phase 
measurements  requires  proper  phase  discrimination,  accurate  determination  of  the  effective  detec¬ 
tion  volume  and  adapt  to  high  concentration  and  glare  near  the  mobile  bed  boundary. 

Technical  Approach 

To  achieve  the  goals  of  this  task,  several  highly  detailed  experiments  have  been  completed  and  are 
being  completed  in  the  controlled  environment  of  a  dust  chamber  in  the  presence  of  a  prototypical 
flow  that  captures  the  essential  features  of  the  rotorcraft  wake  as  it  interacts  with  the  ground.  The 
prototypical  flow  consists  of  a  forced  jet  impinging  on  a  mobile  sediment  bed.  This  configura¬ 
tion  allows  for  the  generation  of  coherent  vortex  rings  embedded  within  a  stagnation  point  flow. 
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Figure  2.16:  a)  Schematic  of  test  setup  for  particle  concentration  calibration,  b)  Measured 
light  sheet  intensity  profile. 

Depending  on  the  forcing  signal,  single  or  multiple  vortex  structures  of  arbitrary  strength  can  be 
generated  and  allowed  to  interact  with  the  stagnating  wall  flow,  providing  for  the  same  essential 
features  of  radial  vortex  stretching  and  vortex/wall  interaction  seen  with  a  rotor  wake.  Quantitative 
measurements  of  the  particle  suspension  and  corresponding  fluid  velocity  are  conducted  through  a 
combination  of  high-speed  imaging  and  two-phase  PIV. 

Results 

Particles  erosion,  suspension  and  turbulence  interaction  with  impinging  forced  jet 

In-situ  concentration  calibration 

Determination  of  the  effective  measurement  volume  is  essential  for  making  quantitative  con¬ 
centration  measurements  of  a  dispersed  phase  when  using  particle-imaging  technique  for  two- 
phase  flows.  Transport  equations  for  both  mass  and  momentum  require  accurate  estimation  of 
particle  volume  fraction  along  with  particle  velocities.  The  effective  measurement  volume,  while 
nominally  determined  by  the  local  light  sheet  thickness,  the  actual  value  depends  also  on  the  dis¬ 
persed  phase  identification  characteristics  used  to  detect  the  particles  (relative  brightness,  size, 
etc.),  and  stray  illumination  such  as  scattering  by  tracer  particles  and  wall  reflections,  necessitat¬ 
ing  use  of  local  calibration  techniques.  A  novel  method  is  being  investigated  where  the  effec¬ 
tive  light  sheet  thickness  is  estimated  using  particle  image  correlation  information  of  free  falling 
dispersed-phase  particles  settling  through  a  tilted  light  sheet  along  with  size  and  brightness  criteria 
of  dispersed  phase  images. 

Moving  glass  particles  through  air  by  fixed  distances  is  challenging,  so  hence  a  different  ap- 
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Figure  2.17:  a)  Particle  brightness  as  a  function  of  position  relative  to  the  light  sheet,  b) 
Measured  particle  size  as  a  function  of  the  position  relative  to  the  light  sheet.  The  nominal 
light  sheet  thickness  as  identified  by  the  beam  profile  camera  is  shown  in  green. 


proach  is  needed  to  estimate  measurement  volume  for  two-phase  PIV  measurements  in  air.  Tests 
have  been  conducted  using  a  high-speed  PIV  imaging  system  working  at  3000  Hz.  The  light  sheet 
is  tilted  with  respect  to  vertical  and  the  camera  is  aligned  perpendicular  to  the  light  sheet  as  shown 
in  Fig.  2.16a.  Glass  particles  of  50-micron  nominal  diameter  are  dropped  into  the  test  section  from 
a  height  of  about  1  m,  so  that  they  enter  the  test  section  moving  at  terminal  velocity.  To  prevent 
external  air  currents  influencing  the  particle  motion,  a  cylindrical  PVC  pipe  is  used  to  shield  the 
falling  particles.  The  terminal  velocity  of  these  particles  is  estimated  to  be  about  1.84  m/s.  The 
estimated  depth  of  focus  is  1.2mm  (105  mm  lens,  f#4)  and  the  field  of  view  is  8  cm2.  The  light 
sheet  thickness  in  the  field  of  view  is  1  mm,  as  shown  in  Fig.  2.16b. 

These  particles  are  imaged  using  a  high-speed  camera  as  they  fall  through  the  light  sheet,  as 
shown  in  the  illustration  in  Fig.  2.16a.  As  particles  move  through  the  light  sheet,  their  scattering 
characteristics  and  hence  their  size  and  brightness  varies.  Shown  in  Fig.  2.17a  is  the  average  peak 
brightness  and  size  of  several  particle  streaks.  When  the  consecutive  images  are  averaged,  the 
particle  paths  appear  as  streaks  in  the  average  image.  The  average  peak  brightness  is  estimated 
from  35  such  particle  streaks  from  the  ensemble,  starting  from  the  first  instance  when  the  particle 
enters  the  light  sheet  until  the  point  where  it  leaves  the  light  sheet.  It  is  evident  that  the  average 
peak  brightness  increases  dramatically  as  the  particles  enter  the  measurement  volume  and  the  av¬ 
erage  size  has  a  more  gradual  increase.  It  can  be  inferred  that  a  particle  image  that  has  its  peak 
brightness  value  above  0.05/o  (where  1Q  is  the  saturation  limit  for  the  camera)  is  within  the  effective 
measurement  volume,  which  in  this  case  turns  out  to  be  2  mm. 

For  a  phase  resolved  low-speed  PIV  system,  it  is  not  possible  to  track  size  and  brightness  of 
individual  particles  as  they  pass  through  the  light  sheet.  In  that  case,  a  PTV  correlation  based 
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Figure  2.18:  a)  Average  correlation  shape  as  a  function  of  the  delay  time  between  images  AT, 
b)  Decay  rate  of  the  correlation  peak  amplitude  as  a  function  of  the  delay  time. 

criteria  has  been  implemented  that  relies  on  decay  of  the  correlation  peak  height  with  increasing 
delay  time  to  identify  the  effective  measurement  volume.  Shown  in  Fig.  4a  is  the  ensemble  av¬ 
erage  correlation  peak  height  of  all  identified  particles  from  5000  images  as  a  function  of  delay 
time.  It  is  evident  that  the  peak  height  drops  because  of  out  of  plane  motion  of  the  particles  with 
increasing  delay  time.  When  these  peaks  are  normalized  with  the  maximum  peak  height  from  the 
lowest  delay  time,  the  decay  trend  can  be  established,  as  shown  in  Fig.  4b.  The  peak  height  drops 
off  dramatically  until  3AT  and  remains  flat  as  the  delay  time  is  further  increased.  This  indicates 
that  these  particles  take  a  time  of  3AT  to  traverse  through  the  measurement  volume  and  knowing 
the  terminal  velocity  of  these  particles  a  priori,  the  effective  light  sheet  thickness  can  be  backed 
out,  which  turns  out  to  be  1.85  mm.  This  is  a  proof  of  concept  study  and  needs  to  be  implemented 
for  current  work  by  matching  the  imaging  conditions  (lighting,  glare  from  the  surface,  FOV  and 
f#)  from  actual  two-phase  experiments. 

Particle  suspension  mechanisms  and  particle-turbulence  interaction 

The  strong  forcing  of  this  flow  produces  a  coherent  (periodic)  structure  that  dominates  the  flow 
behavior.  These  structures  rapidly  erode  the  bed,  generating  specific  bed  forms  that  qualitatively 
alter  the  strength  of  the  suspended  sediment  flux.  The  current  work  has  been  aimed  at  understand¬ 
ing  the  particle-turbulence  coupling  and  net  sediment  transport  at  both  the  early  (denoted  as  ’’flat 
bed”  conditions)  as  well  as  the  late  stages  (denoted  as  ’’eroded  bed”  conditions). 

Given  the  dominance  of  the  forced  periodic  vortex  in  these  flows,  it  is  useful  to  use  ensemble 
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Figure  2.19:  Two  phase  coupling  horizontal  drag  components,  3tzjuDCAu  and  3n/jD(C'Au') 
for  the  flat  bed  case  (a),  and  the  eroded  bed  case  (b).  The  magnitude  is  normalized  by  the  jet 
reference  momentum  flux  through  single  measurement  cell,  pa  UfAyAz. 
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averaging  and  a  triple  decomposition  to  describe  the  flow: 


e  =  e  +  (e)+e'  (2.1) 

where  overbar  represents  the  time-average,  ()  represent  the  coherent  periodic  fluctuations,  and 
represents  the  random  or  stochastic  fluctuations.  The  latter  two  components  by  definition  have  a 
zero  time  average.  Note  that  a  net  ensemble  mean  value  can  be  defined  that  combines  the  effects 
of  the  time-average  and  coherent  fluctuation: 


0  =  0  +  0'  where  0  =  0  +  (0) 


(2.2) 


Once  the  particles  are  in  suspension,  their  interaction  with  the  fluid  can  be  documented  by  estimat¬ 
ing  the  forces  experienced  by  them.  Using  a  point-particle  drag  law  for  dilute  concentrations  of 
particles  and  assuming  a  Stokesian  drag  law,  the  overall  drag  on  particles  in  suspension  is  propor¬ 
tional  to  the  local  particle  concentration  and  the  square  of  slip  velocity: 

D  =  /  3nfjDCAu)  (2.3) 


where 


Au  =  u/  —  Up, 


C  is  the  local  particle  volume  fraction,  D  is  the  particle  diameter  and  /./  is  fluid  dynamic  viscosity. 
Expanding  the  right  hand  side  of  the  drag  equation  using  Reynolds  decomposition  into  mean  and 
fluctuating  components  for  velocity  and  volume  fraction  will  result  in, 


D  =  3nuD 


CAu+(C'Au') 


(2.4) 


The  first  term  in  the  above  equation  is  the  mean  particle-fluid  drag  and  the  second  term  is  the  drag 
fluctuation  due  to  variation  in  slip  velocity.  The  two  contributions  to  the  horizontal  drag  force 
coupling  for  the  flat  and  eroded  bed  conditions  are  highlighted  in  Figure  2.19,  which  are  related  to 
the  turbulent  coupling  between  the  particle  and  fluid  momentum. 

For  the  flat  bed  conditions,  the  suspended  load  interaction  terms  are  confined  quite  close  to 
the  bed  surface  for  both  the  mean  and  fluctuating  contributions  to  the  drag  (Fig.  2.19a),  which  is 
a  signature  of  primary  saltation  behavior  with  little  entrainment  of  particles  into  the  outer  flow. 
There  is  a  small  region  of  positive  drag  on  the  particle  phase  just  beneath  and  slightly  downstream 
of  the  coherent  vortex  structure,  as  particles  are  entering  into  a  high-speed  induced  velocities  by 
the  vortex.  Closer  to  the  bed  and  extended  along  the  majority  of  the  measurement  region,  a  layer  of 
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Figure  2.20:  Suspended  load  transport  for  the  nominally  “flat”  conditions  (top)  and  “eroded 
bed”  conditions  (bottom).  Plots  show  the 

particles  experiencing  a  negative  drag  can  be  observed.  The  fluctuating  component  is  significantly 
larger  in  magnitude  (approximately  20x),  and  also  notably  skewed  towards  positive  contributions. 
This  is  likely  due  to  the  high  degree  of  variability  in  the  position  of  the  vortex  (jitter)  at  these  far 
downstream  positions  caused  by  the  non-linear  interactions  with  the  wall. 

In  the  eroded  bed  case  (Fig.  2.19b),  the  most  notable  difference  is  the  degree  to  which  the 
sediment  is  suspended  away  from  the  wall,  having  significant  interactions  at  distances  up  to  y/R  ~ 
0.3.  The  general  structure  of  the  coupling  is  similar  to  that  observed  for  the  flat  bed  (i.e.  positive 
drag  accelerating  particles  far  from  the  wall,  and  negative  drag  decelerating  particles  near  the  wall), 
with  the  exception  that  the  wave  of  sediment  is  pushed  further  downstream  relative  to  the  vortex 
center. 

The  suspended  load  flux  (i.e.  amount  of  sediment  carried  by  convection  in  suspension)  as  a 
function  of  radial  position  is  calculated  by 

rH 

(<Kr,0)=/  ( C(r,y,t)up(r,y,t))Azdy  (2.5) 

Jo 

and  plotted  in  Figure  2.20  for  the  nominally  flat  and  eroded  bed  conditions.  Although  somewhat 
arbitrary,  here  it  is  presumed  that  any  particle  that  can  be  tracked  by  the  PIV  code  is  sufficiently  far 
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Figure  2.21:  Sediment  bed  profile  evolution  (top),  and  map  of  resulting  denudation  rate 
(bottom). 

from  the  bed  that  they  are  not  undergoing  collisional  reputations  or  minor  saltation  events.  As  we 
will  see  later,  there  is  also  a  significant  bed  load  transport  taking  place.  From  Figure  2.20,  it  can  be 
seen  that  in  both  cases  waves  of  suspended  sediment  are  pushed  in  advance  of  the  passing  vortex. 
For  the  flat  bed  case,  the  suspended  flux  is  continually  increasing  as  a  function  of  radial  distance, 
whereas  for  the  eroded  bed  case,  there  is  little  spatial  variation  in  the  overall  time  average  of  the 
flow. 

The  erosion  of  the  mobile  bed  is  shown  in  detail  for  the  small-particle  low-speed  case  in 
Fig.  2.21.  Tests  were  halted  after  a  maximum  scour  depth  of  10  mm  to  prevent  breakthrough 
of  the  15  mm  thick  bed.  In  general,  the  erosion  starts  just  upstream  of  the  point  of  closest  approach 
of  the  primary  vortex  ( r/R  =  2.0  is  the  point  of  closest  approach  for  the  low-speed  case),  with  a 
crest  and  second  shallower  trough  forming  downstream  of  the  first.  The  wavelength  of  the  bedform 
is  approximately  X/R  -  1.2,  and  varies  little  throughout  the  development  of  the  bed  profile.  A  third 
trough  appears  to  be  forming  beyond  the  second,  but  with  a  greatly  reduced  depth.  The  reduction 
in  depth  with  radial  distance  is  to  be  expected  due  to  both  the  reduced  coherence  of  the  primary 
vortex  as  it  interacts  with  the  ground  surface  (becoming  highly  three-dimensional  and  rapid  dissi¬ 
pation),  as  well  as  due  to  the  deceleration  of  the  wall  jet  due  to  the  axisymmetric  expanding  flow. 
Once  the  initial  trough  begins  to  form,  scouring  occurs  progressively,  and  moves  the  troughs  and 
crests  downstream  while  keeping  the  lee  slope  profile  into  the  first  trough  approximately  constant. 
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The  denudation  rate  is  defined  simply  as  the  rate  of  removal  of  the  bed  surface,  and  is  calculated 
locally  by  subtracting  the  difference  in  elevation  of  the  bed  at  successive  times  (see  Fig.  2.21).  This 
clearly  reveals  two  maxima  in  surface  erosion  rate  located  about  the  center  of  the  nascent  forming 
troughs  at  the  earliest  times  ( r/R  ~  2.3  and  r/R  ~  3.5).  Although  it  seems  clear  that  the  forcing 
by  the  primary  vortex  is  responsible  for  the  initiation  of  the  bed  forms,  due  to  the  correlation  of 
the  periodic  Reynolds  stress  with  the  initial  erosion  maxima,  it  is  not  yet  clear  from  the  evidence 
we  have  as  to  what  causes  the  second  erosion  maxima.  It  may  possibly  that  this  results  from  an 
increase  in  the  stochastic  component  of  the  turbulent  stresses  caused  by  the  breakdown  and  three- 
dimensionalization  of  the  primary  vortex,  which  occurs  in  this  region.  Soon  after  formation  of  the 
bedforms,  the  peak  erosion  rates  shift  to  the  stoss  slope  with  similar  erosion  rates  on  both  dunes,  in 
contrast  to  the  marked  asymmetry  shown  between  the  two  dunes  in  the  initial  erosion  conditions. 
Minimal  erosion  (or  even  deposition)  is  evident  on  the  lee  slopes. 

Finally,  the  net  erosion  rate  over  the  measured  field  of  view  is  determined  by  total  rate  of  area 
removed,  as  plotted  in  Fig.  2.22.  It  can  be  seen  that  the  early  times  ( t  <  30  seconds,  correspond¬ 
ing  to  the  flat  bed  case),  the  net  erosion  rate  is  comparable  to  or  greater  than  what  was  observed 
under  eroded  bed  conditions  ( t  «  80  second).  This  interesting  as  the  suspended  flux  in  the  eroded 
bed  conditions  shown  in  Fig.  2.20  carries  approximately  10  times  the  amount  of  sediment  in  com¬ 
parison  to  the  flat  bed  conditions.  This  is  indicative  of  the  flat  bed  being  dominated  by  bed  load 
transport,  and  the  role  that  the  bed  forms  play  in  ejecting  sediment  into  the  upper  regions  of  the 
flow.  The  second  notable  point  in  Fig.  2.22b  is  the  fact  that  the  sediment  transport  scales  with  the 
number  of  vortex  structures  that  have  occurred.  Once  scaled,  the  timing  of  the  events  in  the  erosion 
curves  similar  transient  behavior,  underscoring  the  fact  that  the  erosion  is  being  dominated  by  the 
vortex  events. 
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Figure  2.22:  a)  Dimensional  erosion  rates,  b)  Erosion  depth  per  vortex  impact. 
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Task  2.4 


Large-Eddy  Simulation  of  the  Interaction  Between  Wake  Vortices  and  Ground 

Investigator(s):  Ugo  Piomelli 

Institution/Department:  Queens  University,  Dept,  of  Mechanical  and  Materials  Engineering 
Contact  email(s):  ugo@me.queensu.ca 

This  task  was  part  of  the  original  MURI  proposal,  and  similar  work  is  currently  being  per¬ 
formed  under  separate  AFOSR  funding  by  Professor  Piomelli  at  Queens  University  in  Canada. 
Further  information  on  this  work  can  be  obtained  directly  from  Professor  Piomelli. 
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Task  2.5 


Two-Phase  Large  Eddy  Simulation  Using  the  Mesoscopic  Eulerian  Formal¬ 
ism 

Investigator(s):  Kyle  Squires 

Institution/Department:  Arizona  State  University,  Mechanical  and  Aerospace  Engineering 
Graduate  Student(s):  Dennis  Dunn 
Contact  email(s):  squires@asu.edu 

Background  and  Technical  Challenges 

The  brownout  problem  encompasses  multiple  scales,  from  the  largest,  comparable  to  the  dimen¬ 
sions  of  the  rotorcraft,  to  the  smallest,  comparable  to  the  sediment  particle  size.  This  extended 
range  leads  to  complex  effects  that  influence  the  processes  of  entrainment,  deposition,  and  resus¬ 
pension  of  the  dispersed  sediment.  Once  suspended,  sediment  particle  interactions  occur  with  the 
coherent  wake  vortices  characterizing  the  rotor  flow,  and  with  the  finer  scale  turbulence  generated 
near  the  ground  as  it  interacts  with  the  sediment  bed.  Factors  that  influence  transport  charac¬ 
teristics  of  the  resulting  two-phase  flow  include  particle-turbulence  interactions,  particle-particle 
interactions,  particle  interactions  with  the  bed  (e.g.,  entrainment,  suspension,  deposition),  and  a 
momentum  coupling  that  modifies  the  carrier  flow. 

In  general,  the  primary  technical  challenges  can  be  sub-divided  into  those  residing  at  the  mi¬ 
croscale  (essentially  the  scale  of  a  single  particle)  and  the  mesoscale  (essentially  the  scale  com¬ 
prised  of  a  large  ensemble  of  particles  with  dimensions,  for  example,  comparable  to  the  charac¬ 
teristic  boundary  layer  thickness  of  the  two-phase  flow).  Specific  technical  challenges  include 
simulation  or  measurement  of  the  detailed  processes  of  particle  lift-off  from  a  sediment  bed,  the 
role  of  turbulence  influencing  dispersed-phase  transport,  particle-particle  interactions,  and  the  in¬ 
fluence  of  particle  momentum  exchange  on  properties  of  the  carrier  phase  turbulent  flow. 

Modeling  the  transport  of  a  second  phase  of  small,  dispersed  particles  within  a  carrier-gas- 
phase  fluid  flow  has  been  successfully  simulated  in  this  MURI  task  by  using  Eulerian-Lagrangian 
approaches  in  which  Discrete  Particle  Simulation  (DPS)  is  used  to  compute  properties  of  the  parti¬ 
cle  phase.  Eulerian-Lagrangian  approaches  have  been,  and  will  remain,  an  invaluable  tool  though 
are  often  constrained  to  flows  with  very  small  volume  fractions  and  low  mass  loadings.  Addition¬ 
ally,  Lagrangian  particle  tracking  approaches  become  more  demanding  in  terms  of  inter-processor 
communication  and  load  balancing  for  non-uniformly  distributed  particle  concentrations  and  can 
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be  very  challenging  to  model  large  particle  ensembles.  Eulerian-based  methods  overcome  these 
challenges  and  present  a  more  straightforward  implementation  on  a  solution  space  that  is  divided 
among  parallel  CPU’s. 

The  above  considerations  provide  the  main  motivation  for  the  development  of  “two-fluid” 
Eulerian-Eulerian  methods  in  which  the  dispersed  phase  is  modeled  as  a  second  fluid.  These 
approaches  provide  a  sound  framework  for  modeling  flows  with  large  particle  concentrations  and 
a  range  of  complex  interactions,  e.g.,  particle  collisions,  breakup,  coalescence,  etc.  Eulerian-based 
approaches  also  present  computational  simplifications  in  that  the  same  basic  framework  for  so¬ 
lution  of  the  field-based  equations  can  be  applied  as  used  to  model  fluid  motion  and  overcome 
inter-processor  communication  and  load  balancing  problems  that  can  arise  for  non-uniformly  dis¬ 
tributed  particle  concentrations.  Finally,  Eulerian-based  approaches  provide  a  pathway  to  true, 
two-fluid  Large  Eddy  Simulation  of  multiphase  flows.  The  approach  that  has  been  implemented 
and  tested  for  Eulerian-based  predictions  of  the  rotor-vortex  boundary  layer  is  based  on  Condi¬ 
tional  Quadrature-Based  Method  of  Moments  (CQMOM).  This  latter  method  is  being  used  to 
model  dispersed  phase  transport  and  particle  interactions  and  is  being  assessed  against  Lagrangian- 
based  modeling  strategies  that  have  been  successfully  implemented  in  earlier  phases  of  this  task. 

Technical  Approach 

The  starting  point  for  Eulerian-based  approaches  is  introducing  a  probability  density  function 
(PDF),  /,  that  describes  the  probability  of  finding  a  particle  at  time  t  and  position  x  possess¬ 
ing  velocity  v,  diameter  d,  temperature  T,  etc.  These  properties  are  parameters  of  the  PDF, 
f(t,x,v,d, T, ...),  where  t  and x  are  independent  variables  and  the  remaining  variables  {v,d, T,...} 
define  internal  state  space  coordinates.  This  approach,  in  turn,  corresponds  in  practical  compu¬ 
tations  to  cell-localized  quantities  representing  the  likelihood  of  a  specific  coordinate  set  (e.g., 
(vo,do,  7o, ...})  existing  in  that  cell.  For  finite-inertia  particles  whose  resistance  to  the  external 
flow  is  characterized  by  the  Stokes  number,  St,  the  individual  particles  in  the  ensemble  may  cross 
paths  and  exhibit  diverse  velocities.  Consequently,  the  values  of  the  internal  variable  set  (e.g., 
velocity,  diameter,  and  temperature)  at  a  given  point  in  space  can  assume  more  than  one  value, 
again  owing  to  the  fact  that  particles  possessing  different  velocities,  diameters,  or  temperatures 
can  occupy  the  same  cell  over  which  averages  are  computed.  It  is  thus  important  for  modeling 
particle-laden  flows  to  account  for  multi-valued  solutions.  This  goal  is  accomplished  in  the  con¬ 
text  of  the  PDF  by  introducing  the  transport  equation  for  /  that  for  a  mono-disperse,  isothermal 
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flow  ( d  and  T  are  constant)  is  a  specific  version  of  the  Boltzmann  Equation,  i.e., 

3/  9/  d 

¥  ¥  ¥'(v/)_c 


(2.6) 


where  C  is  a  collision  source  term  and  v  is  particle  acceleration  (dictated  by  the  drag  force  in 
work  performed  to  date).  The  above  approach  is  general  in  that  the  Navier-Stokes  equations  for 
continuous  fluids  are  also  derived  from  Eq.  2.6.  In  cases  where  particle  trajectories  must  be  allowed 
to  cross,  as  occurs  in  flows  with  finite-inertia  particles  that  do  not  follow  the  fluid  streamlines,  the 
assumption  that  at  a  given  location  in  space  the  particle  velocity  v(t,x)  is  single-valued  is  not  made. 

For  practical  simulation  strategies,  it  is  not  feasible  to  directly  solve  Eq.  2.6  and  instead  alterna¬ 
tive  equations  are  derived  by  taking  moments  of  Eq.  2.6,  leading  to  transport  equations  of  the  same 
form  as  those  obtained  for  continuous  fluids  governing  conservation  of  mass,  momentum,  energy, 
etc.  Without  loss  of  generality,  the  moments  of  /  may  be  expressed  in  the  following  manner 


Mlj,k(^x)  =  f  V\V2V3  f(t,x,v)d\ 


(2.7) 


where  i,j,k  are  the  chosen  integer  exponents  used  to  skew  the  PDF  integration  by  the  internal 
particle  velocity  variables  v  =  {vi ,  V2 ,  V3  },  and  y  =  i  +  j  +  k  is  the  sum  of  the  exponents  that  simply 
indicates  the  order  of  the  moment.  This  skewed  integration  of  the  PDF  relies  on  carefully  chosen 
integer  values  (>  0)  for  the  exponents  i,j,k  that  produce  the  optimal  set  of  moments.  Applying 
this  moment  transformation  to  the  PDF  transport  in  Eq.  2.6,  produces  a  general  set  of  integro- 
differential  equations,  i.e.. 


|M^+  Vx-Afl+l  =  z*  v-AfY-1  +  C 
Unsteady  Convective  Acceleration  Source 


(2.8) 


where  the  operation  i  *  v  =  {iv\ , ./v’2 ,  A'V’3 } .  For  clarity,  the  equation  set  described  by  Eq.  2.8  shows 
only  the  order  of  the  moment  y  with  the  subscripts  removed.  The  acceleration  term  is  a  source  term 
that  arises  from  drag,  gravity,  or  other  body  forces  and  involves  moments  of  order  y—  1,  while  the 
convective  term  always  involves  moments  of  order  y+  1.  Equations  2.6  through  2.8  are  exact 
representations,  but  to  close  the  higher  order  convective  moments  (y+  1).  A  Gaussian  quadrature 
approximation  is  applied  to  Eq.  2.7,  i.e., 

P 

~  £  UlaVlW^a{x,t)  Y  =  0, 1,2,...  (2.9) 

0C=1 

where  |3  is  the  total  number  of  discrete  nodes  allocated  to  approximate  the  function  /  as  a  sum  of 
Dirac  Delta  functions.  A  single  “node”  describes  a  complete  set  of  abscissas  and  associated  weight: 
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{U,  V,  W, (|)}.  The  D  =  3  variables  U,  V,  W  are  the  three-dimensional  particle  velocity  abscissas  and 
<|)  the  associated  weight  for  the  a’th  quadrature  node.  In  total,  this  requires  K  =  (3D  •  (D  +  1)  total 
moment  equations  to  provide  closure  for  the  set  of  abscissas  and  weights  in  CQMOM. 

By  this  method,  /  is  reduced  to  a  (3-valued  model  in  which  up  to  (3  distinct  instances  of  the 
nodes  may  be  present  in  any  cell.  The  Gaussian  quadrature  approximates  the  area  under  /  such 
that  the  moment  set  retains  information  for  a  multi-valued  solution.  Numerically  solving  these 
moment  equations  is  challenging,  but  can  also  be  extremely  effective  and  useful,  even  using  the 
approximation  of  (3  =  2  employed  in  the  present  model.  It  should  be  noted  that  for  Eulerian  models 
that  employ  a  single-valued  assumption,  such  as  the  Navier-Stokes  equations,  only  one  realization 
of  each  internal  velocity  variable  ((3=1)  would  be  stored  at  every  grid  point,  consistent  with  the 
fact  that  the  ability  to  account  for  particle  trajectory  crossing  is  not  required. 

Together,  the  moment  set  M  and  nodes  {U,  V,  W,  <|)}  transport  the  necessary  quantities  for  mod¬ 
eling  the  flow  while  offering  the  ability  to  impose  boundary  conditions  and/or  body  forces  in 
a  physically  feasible  manner.  This  particular  version  of  the  CQMOM  model  was  proposed  by 
Yuan  and  Fox*  and  is  employed  in  the  present  effort  because  it  offers  the  ability  to  represent  a 
three-dimensional  multi-valued  velocity  field  for  a  dispersed  Eulerian  phase,  and  it  has  yielded 
promising  results  for  monodisperse  particle  phase  distributions  in  turbulent  flows. 

Results 

The  CQMOM  model  has  been  incorporated  into  the  DNS  fluid  solver  that  has  been  employed 
for  Lagrangian/DPS  studies  previously  during  this  MURI  task.  In  total,  two  phases  are  being 
modeled:  gas-phase  carrier  flow,  and  the  dispersed  particle  phase  is  being  modeled  in  two  different 
ways.  Lagrangian  particles  are  modeled  using  DPS,  and  an  Eulerian  treatment  of  the  particle  phase 
is  modeled  using  CQMOM  (referred  to  as  “Eulerian  particles”  below).  For  comparative  purposes, 
the  two  particle  phases  run  concurrently  and  do  not  interact  with  one  another,  and  with  one-way 
coupling  to  the  fluid  that  in  turn  ensures  each  particle  treatment  experiences  the  same  forcing  from 
the  fluid.  The  Stokes  number  and  inter-collision  states  (enabled/disabled)  are  identical  between  the 
Eulerian  and  Lagrangian  particle  simulations  to  enable  side-by-side  comparisons.  The  following 
figures  may  compare  all  three  models  as  appropriate.  Additionally,  instantaneous  snapshots  and 
time-averaged  statistics  are  available  for  selected  2-D  planes  within  the  computational  domain. 

*Fox,  R.  O.  “Higher-order  Quadrature-based  Moment  Methods  for  Kinetic  Equations,”  Journal  of  Com¬ 
putational  Physics,  228,  No.  20,  2009,  pp.  7771-7791;  Yuan,  C.  and  Fox,  R.  O.,  “Conditional  Quadrature 
Method  of  Moments  for  Kinetic  Equations,”  Journal  of  Computational  Physics,  Vol.  230,  No.  22,  201 1,  pp. 
8216-8246. 
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Figure  2.23:  Instantaneous  streamwise  velocity  at  half-span,  St= 4,  without  inter-particle  col¬ 
lisions. 


The  spanwise  z  direction  of  the  computational  domain  is  periodic,  while  the  upper  y  boundary 
employs  a  Neumann  boundary  condition,  and  the  lower  y  boundary  acts  as  a  reflective  wall.  The 
streamwise  inlet  x  direction  maintains  a  fixed  weight  distribution  and  the  outlet  x  direction  also 
uses  Neumann  boundary  conditions. 

Two  separate  cases  are  considered.  The  first  uses  the  turbulent  boundary  layer  (TBL)  carrier- 
flow  field  (without  coherent  vortices  that  impinge  into  the  boundary  layer)  to  drive  the  lighter 
particles  (St  =  4)  while  neglecting  inter-particle  collisions.  Figure  2.23  shows  an  instantaneous 
snapshot  over  about  25%  of  the  extent  of  the  streamwise  domain  for  the  fluid  flow  boundary  layer 
and  both  models  of  the  particle  phases,  both  at  the  same  instant  in  time.  The  shear  from  the  no¬ 
slip  condition  on  the  wall  creates  complex  vortical  structures  in  the  boundary  layer  even  without 
coherent  vortices  impinging  onto  the  flow.  The  results  in  the  figure  shows  that  the  lighter  particles 
follow  the  fluid  flow  field  fairly  closely  with  the  largest  velocity  discrepancies  occurring  in  regions 
of  rapid  change  (large  velocity  gradients)  or  small  localized  perturbations,  owing  to  the  inertia- 
induced  particle  lag.  Importantly,  the  figure  also  shows  good,  qualitative  agreement  between  the 
Lagrangian  and  Eulerian  models  of  the  particle  motion,  in  turn  providing  a  key  validation  of  the 
Eulerian-based  approach  for  predicting  the  dispersed  phase  properties. 

Unlike  the  streamwise  velocity,  both  the  wall-normal  and  spanwise  velocities  (v2  and  V3,  re¬ 
spectively)  possess  zero  (for  V3)  or  nearly  zero  (for  V2)  long  time-averages  and  deviations  from  the 
fluid  flow  will  be  more  apparent,  as  shown  in  Fig.  2.24.  Both  the  Lagrangian  and  Eulerian  particle 
treatments  exhibit  large  regions  of  correlated  motion  that  is  driven  by  the  interaction  of  the  parti¬ 
cles  with  the  motion  of  the  fluid  flow.  The  particle  velocity  magnitudes  shown  in  Figs.  2.24b  and 
2.24c  are  about  half  as  large  as  those  in  the  fluid  shown  in  Fig.  2.24a,  though  again  exhibit  very 
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A 


(a)  Carrier-fluid  (b)  Eulerian  particles 


(c)  Lagrangian  particles 


Figure  2.24:  Instantaneous  wall-normal  velocity  in  a  wall-parallel  plane  at  y  «  0.31,  St  =  4, 
without  inter-particle  collisions. 


(a)  Contours  of  wall-normal  fluid  velocity,  indicating  locations  of  the  velocity  profiles. 
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(b)  Profiles  of  the  wall-normal  V2  mean  velocity. 


Figure  2.25:  Time-averaged  wall-normal  velocity  at  half-span,  St= 4,  without  inter-particle 
collisions. 


similar  structure.  Although  the  particles,  in  general,  follow  the  fluid,  the  important  validation  of 
CQMOM  is  realized  because  the  two  particle  models  are  in  good  qualitative  agreement  with  one 
another. 

The  final  figure  for  the  TBL  case  in  Fig.  2.25b  shows  time-averaged  wall-normal  velocities 
using  line  plots  for  comparison.  The  plane  lies  at  the  half-span  and  has  accumulated  statistics  that 
enable  some  assessment  of  the  two  models.  The  results  show  that  CQMOM  generally  predicts 
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(b)  Profiles  of  the  wall-normal  V2  mean  velocity. 


Figure  2.26:  Time-averaged  wall-normal  velocity  at  half-span,  St  =  16,  with  inter-particle 
collisions. 


larger  magnitudes  than  DPS,  though  this  outcome  is  an  artifact  of  the  relatively  short  statistical 
sample  currently  available.  In  sparsely-populated  regions  where  cells  are  frequently  devoid  of 
Lagrangian  particles,  the  DPS  solution  will  yield  lower  averages,  which  also  contributes  to  the 
much  nosier  profiles  shown  in  the  figure.  In  denser  regions  there  is  better  agreement  between  the 
CQMOM  and  DPS  models.  Additionally,  the  Eulerian  particle  velocities  are  generally  different  in 
magnitude  than  the  fluid  owing  to  their  inertia. 

The  second  case  with  larger-inertia  particles,  St  =  16,  differs  from  the  first  in  many  aspects:  a 
time-periodic  vortex  impinges  into  the  boundary  layer  near  x  «  20  and  with  a  larger  Stokes  number 
that  results  in  greater  particle  trajectory  crossing,  and  an  inter-particle  collision  model  is  employed. 
For  these  larger  inertia  particles  the  fluid  and  particle  velocities  can  deviate  significantly.  Similar 
time- averaged  statistical  measures  are  shown  in  Fig.  2.26b,  sampled  over  3.85  time  units,  and  show 
general  qualitative  agreement  between  CQMOM  and  DPS,  though  with  the  larger  Stokes  number 
the  particles  do  not  always  conform  as  closely  to  the  fluid  velocities. 

The  Lagrangian  particles  are  initially  randomly  seeded  throughout  the  domain,  which  impacts 
the  current  averages  sampled  over  a  relatively  short  period  because  the  cells  devoid  of  particles 
artificially  lower  the  average.  For  instantaneous  velocities  and  densities,  the  DPS  model  some¬ 
times  has  difficulty  describing  the  larger-scale  structures  where  an  insufficient  quantity  of  particles 
results  in  empty  cells  that  are  given  zero  velocity  in  the  initialization.  This  issue  can  cause  La- 
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(a)  Eulerian  density  M°  (b)  Lagrangian  particles  per  cell. 


Figure  2.27:  Instantaneous  particle  density  in  a  wall-parallel  plane  (y  ~  0.31),  57=16,  includ¬ 
ing  inter-particle  collisions. 


(c)  Lagrangian  particle  field 


Figure  2.28:  Time-averaged  wall-normal  velocity  V2  in  a  wall-parallel  plane  (y  «  0.31),  57=16, 
including  inter-particle  collisions. 

grangian  time-averages  to  require  more  samples  than  may  be  necessary  than  the  Eulerian  counter¬ 
part.  The  CQMOM  model  tracks  a  density  quantity  M°  that  was  seeded  more  uniformly  no  larger 
than  an  initial  value  of  0.3  for  ease  of  implementation  and  an  initial  assessment  of  the  model. 

The  dispersion  and  accumulation  of  particles  through  interaction  with  the  turbulent  carrier 
flow  quickly  leads  to  sparse  regions  with  M°  0  and  higher  concentration  regions  where  M°  «  1 
or  more,  showing  that  the  concentrations  may  achieve  at  least  three  times  their  initial  values.  Com¬ 
parisons  can  be  made  between  the  density  fields  of  the  two  particle  models  (the  incompressible 
fluid  field  is  omitted),  as  shown  in  Fig.  2.27.  The  Lagrangian  density,  on  the  other  hand,  is  mea¬ 
sured  in  “particles  per  bin”  and  yields  very  similar  structures  for  the  particle  densities,  which  are 
discernible  upon  careful  inspection  of  Figs.  2.27a  and  2.27b.  Selected  regions  are  highlighted  to 
showcase  the  similarity  of  the  structures,  which  provides  an  important  and  powerful  validation  of 
the  CQMOM  model. 
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(c)  Lagrangian  particle  field 

Figure  2.29:  Time-averaged  spanwise  velocity  V3  in  wall-parallel  plane  (y  «  0.0232),  St  =  16, 
with  inter-particle  collisions. 
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Figure  2.30:  Instantaneous  particle  density  in  a  wall-parallel  plane  (y  ~  0.0232),  St  =  16,  with 
inter-particle  collisions. 

Good  agreement  is  observed  between  the  two  particle  models  in  most  cases.  Small  discrep¬ 
ancies  may  appear  in  the  results  showing  instantaneous  quantities  where  the  DPS  method  lacks 
sufficiently  large  sample  sizes  in  some  bins,  giving  statistical  outliers  more  weight.  CQMOM  is 
essentially  showing  a  statistical  representation  of  an  arbitrarily  large  quantity  of  pseudo-particles 
averages  as  a  continuous  medium.  The  wall-normal  particle  velocities  are  very  similar,  show¬ 
ing  similar  structures  in  Fig.  2.28b,  even  with  the  sparse  Lagrangian  particle  concentrations  in 
Fig.  2.28c,  showing  smaller  values  overall. 

Figure  2.29  similarly  shows  the  spanwise  velocity  field  near  the  wall,  and  shows  that  CQMOM 
predicts  slightly  greater  velocity  magnitudes  than  occur  in  the  fluid,  though  maintains  good  qual¬ 
itative  agreement  with  the  structural  features  in  the  flow  between  Figs.  2.29b  and  2.29c.  For  both 
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the  wall-normal  and  spanwise  velocities  of  Figs.  2.28  and  2.29,  respectively,  CQMOM  predictions 
in  (b)  show  larger  magnitudes  but  also  show  structures  more  correlated  to  that  of  the  Lagrangian 
particles  (c)  than  to  that  of  the  fluid  flow  (a). 

Figure  2.30  shows  a  view  of  the  instantaneous  particle  density  over  the  entire  x-z  domain  near 
the  wall  (y  ~  0.0232)  where  the  initial  seeding  was  M°  ~  0.15.  These  results  were  previously  not 
realizable  using  DPS  due  to  an  insufficient  number  of  particles  in  that  region.  The  near  wall  plane 
shown  in  Fig.  2.30  shows  interesting  “horseshoe”  structure  formations  in  the  density  field. 
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Task  2.6 


Effects  of  Model  Scaling  on  Sediment  Transport 

Investigator(s):  Anya  Jones 

Institution/Department:  University  of  Maryland/Department  of  Aerospace  Engineering 
Graduate  Student(s):  Gino  Perrotta,  Mark  Glucksman-Glaser 
Contact  email(s):  arjones@umd.edu 

Background  and  Technical  Challenges 

Several  experiments  have  been  conducted  to  study  rotor  aerodynamics  in  ground  effect  and  the  par¬ 
ticle  phase  transport  dynamics  that  underly  the  problem  of  rotorcraft  brownout  (e.g.,  Tasks  1.1  and 
2.3).  Because  full-scale  rotor  experiments  are  impractical,  the  research  has  relied  on  laboratory- 
scale  experiments  in  a  controlled  environment  to  investigate  the  various  phenomena.  However, 
although  the  bulk  of  the  ongoing  work  currently  relies  on  small  scale  experiments,  the  effects  of 
scaling  on  the  aerodynamics  (carrier  phase)  and  particle  transport  (dispersed  phase)  in  a  complex 
and  tightly-coupled  dual-phase  flow  such  as  brownout  are  not  yet  understood. 

The  brownout  problem  requires  that  both  the  aerodynamics  and  the  sediment  transport  be 
scaled  together,  and  so  the  many  dimensionless  similarity  parameters  associated  with  both  the 
aerodynamics  and  the  sediment  must  be  matched  simultaneously  to  obtain  a  truly  scaled  dual¬ 
phase  flow.  Because  this  requirement  is  nearly  impossible  to  meet,  a  sensitivity  study  can  be  used 
to  gain  confidence  that  laboratory-scale  experimental  results  are  applicable  to  the  full-scale  rotor- 
craft  brownout  problem.  The  focus  of  this  task  is  to  identify  aerodynamic  and  sediment  transport 
similarity  parameters  that  can  be  used  to  relate  laboratory-scale  studies  of  rotor-induced  sediment 
mobilization  and  uplift  to  full  scale  brownout  events.  The  classical  aerodynamic  and  sediment 
transport  similarity  parameters  have  been  evaluated,  and  the  sensitivity  of  the  dual-phase  flow  ob¬ 
served  in  the  laboratory  to  these  parameters  has  been  assessed.  The  current  experiments  are  being 
performed  in  water  rather  than  air  to  allow  for  a  broad  parameter  space  for  sensitivity  studies  by 
increasing  the  viscosity  of  the  fluid  and  so  changing  the  the  time  and  length  scales  of  the  flow. 

The  primary  objectives  of  this  task  are:  1.  To  determine  which,  if  any,  of  the  classical  aerody¬ 
namic  and  sediment  transport  similarity  parameters  best  characterize  the  scaling  of  brownout;  2. 
To  identify  new  scaling  parameters;  3.  To  evaluate  the  sensitivity  of  brownout  (severity,  sediment 
uplift/transport  mechanisms,  sediment  dispersion/suspension,  etc.)  to  scaling  over  a  wide  range  of 
model  scales. 
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Technical  Approach 

Experiments  have  been  performed  on  a  single-bladed  rotor  (thrust-matched  to  previous  experi¬ 
ments  in  air  performed  under  Task  2.2)  operating  in  ground  effect  above  a  sediment  bed  in  a  water 
tank  of  dimensions  1.2  x  1.2  x  1.2  m.  Several  types  of  sediment  particles  were  used  to  provide 
a  wide  range  of  scaling  parameter  values.  The  experiments  have  made  use  of  particle  image  ve- 
locimetry  (PIV)  to  quantify  the  carrier  flow  as  well  as  the  quantity,  location  and  trajectory  of  the 
sediment  particles.  The  specific  goals  were:  1.  Quantify  the  behavior  the  carrier  phase  using  time- 
resolved  PIV  measurements  of  the  rotor  wake  to  evaluate  the  relationship  between  unsteady  flow 
phenomena  (e.g.,  a  vortex  passing  over  the  sediment  bed)  and  sediment  uplift  and/or  transport.  2. 
Assess  the  relationship  between  sediment  uplift  in  a  scaled  brownout  flow  environment  to  classical 
aerodynamic  and  sediment  transport  scaling  parameters.  3.  Investigate  and  evaluate  new  potential 
scaling  parameters  for  the  dual-phase  flow  produced  by  a  rotor  that  is  operating  in  ground  effect 
over  a  mobile  sediment  bed. 

Results 

The  wake  generated  by  a  rotor  in  ground  effect  is  highly  unsteady  and  cannot  be  adequately  char¬ 
acterized  near  the  ground  by  a  time-averaged  flow  field.  Figure  31(a)  shows  a  time-averaged 
velocity  field  in  the  region  where  the  rotor  wake  impinges  on  the  ground.  Figure  31(b)  shows  an 
instantaneous  velocity  field  of  the  same  region  in  which  cross-sectional  elements  of  the  helical  tip 
vortex  can  be  seen  impinging  on  the  ground  plane  and  producing  large  velocity  transients.  Previ¬ 
ous  experiments  in  air  (Task  2.2)  have  suggested  that  the  unsteady  local  disturbances  caused  by 
the  passage  of  these  vortices  initiates  the  sediment  uplift  cycle.  This  latter  finding  has  been  sup¬ 
ported  by  flow  visualization  experiments  in  water  that  have  been  conducted  under  the  current  task. 
Sediment  uplift  has  been  observed  primarily  just  below  and  ahead  of  a  vortex. 

One  research  challenge  is  that  of  characterizing  the  unsteady  features  of  the  rotor  wake.  To 
enable  comparisons  between  different  cases  and  different  studies,  the  wake  characteristics  were 
summarized  in  a  concise  set  of  scalar  quantities,  the  primary  flow  feature  of  which  is  the  rotor  tip 
vortex,  which  can  be  described  by  three  quantities:  the  total  circulation  contained  in  the  rotating 
flow,  Tv,  the  maximum  tangential  flow  velocity  around  the  vortex  center,  Vemax,  and  the  vortex 
core  radius,  rc.  Given  a  model  for  the  variation  in  tangential  velocity  along  the  vortex  radius,  any 
one  of  the  three  parameters  can  be  found  in  terms  of  the  other  two. 

A  vortex  tracking  program  was  developed  and  used  to  identify  the  locations  of  the  vortices  in 
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x  rotor  radius  x  rotor  radius 

(a)  Time-averaged  flow  field  (20  revolutions)  (b)  Instantaneous  flow  field 

Figure  2.31:  Contours  of  magnitude  of  total  velocity  and  velocity  vectors  for  a  rotor  wake 
impinging  on  a  ground  plane.  85  mm  radius  rotor  hovering  at  300  rpm.  Rotor  is  one  radius 
above  ground  and  the  field  of  view  is  approximately  1-2  radii  from  axis  of  rotation  and 
approximately  1  radius  from  the  ground  plane. 

each  flow  field.  Each  of  the  identified  vortices  in  each  image  was  then  measured  in  terms  of  VQmax, 
rv,  and  rc.  This  was  accomplished  by  fitting  the  Lamb-Oseen  vortex  profile  to  a  scatter  plot  of  the 
local  tangential  velocities  compared  to  radial  distance  from  the  vortex  center.  The  best  fit  of  this 
model  to  the  data  provided  a  repeatable  quantification  of  the  size  and  strength  of  each  vortex  flow, 
which  in  turn  allows  for  quantification  of  the  flow  field  in  terms  of  unsteady  characteristics. 

When  compared  to  the  rotor  speed,  two  of  the  mean  vortex  parameters,  Eemax  and  Tv,  form 
nearly  linear  trends  (see  Figure  2.32).  This  observation  agrees  with  intuition  because  the  energy 
in  the  flow  increases  with  the  rotor  speed,  and  the  increased  energy  would  create  higher  rotational 
speed  and  higher  circulation  in  the  vortices.  The  highly  non-linear  trend  for  rc  warrants  more 
discussion.  The  core  radius  was  found  to  decrease  as  the  rotor  speed  is  increased  from  the  lowest 
speed  cases  to  the  middle  of  the  range.  After  that,  the  radius  increased  with  the  rotor  speed  out  to 
the  highest  speed. 

Two  primary  phenomenon  are  in  conflict  here.  As  the  rotor  speed  increases,  the  lift  of  the 
blade  is  also  increased.  High  lift  creates  a  stronger  tip  vortex,  which  more  closely  resembles  an 
irrotational  vortex.  At  the  extreme  limit  of  this,  the  peak  velocity  would  occur  at  the  vortex  center. 
Approaching  this  limit  means  that  the  peak  swirl  velocity  occurs  closer  to  the  center,  and  so  the 
stronger  vortices  have  a  smaller  core  radius.  At  the  same  time,  higher  speed  flow  experiences 
increased  shear  stress.  This  stress  can  promote  turbulence  that  will  increase  the  rate  of  diffusion, 
spreading  the  vorticity  to  the  surrounding  flow  and  increasing  the  measured  core  radius.  The  trend 
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Figure  2.32:  Wake  characteristics  plotted  with  respect  to  rotor  speed 
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observed  for  rc  suggests  that  the  core-reducing  effect  is  dominant  at  low  rotor  speeds  and  the 
core-expanding  effect  is  dominant  at  high  rotor  speeds. 

Given  time-resolved  images  of  the  two  phase  flow,  the  “severity”  of  the  particle  field  in  the  dis¬ 
persed  phase  was  quantified  by  counting  the  total  number  of  sediment  particles  suspended  in  the 
field  of  view.  This  quantity  provides  a  simple  measure  of  the  intensity  of  the  brownout,  but  could 
be  further  modified  to  also  include  the  height  and/or  velocity  of  the  particles  above  the  bed.  A  min¬ 
imum  threshold  and  median  filter  were  used  to  remove  most  of  the  nylon  tracer  particles  from  the 
images.  Using  Khalitov’s  search  function,  Gaussian-like  peaks  in  image  intensity  were  identified 
as  individual  particles,  and  another  minimum  threshold  was  then  used  to  eliminate  any  data  that 
were  not  part  of  the  Gaussian  peaks.  A  sediment  particle  was  defined  to  have  a  certain  minimum 
size  and  minimum  mean  brightness  in  the  image;  data  below  threshold  values  were  eliminated. 
A  maximum  particle  size  was  implemented  to  eliminate  groups  of  sediment  particles  that  were 
not  uplifted  but  were  identified  by  the  Gaussian  search.  The  final  sediment  search  function  was 
controlled  by  the  following  five  parameters:  median  filter  width,  minimum  threshold,  minimum 
and  maximum  sediment  size,  and  minimum  sediment  brightness.  For  two  images,  each  visible 
sediment  particle  was  manually  identified.  Values  for  each  search  parameter  were  tested  using  a 
Monte  Carlo  method,  where  results  were  scored  for  accuracy  compared  to  the  manually  identified 
particle  locations.  This  process  was  run  separately  for  the  steel  and  glass  sediment.  Given  a  time 
series  of  PIV  images,  a  record  of  sediment  uplift  as  a  function  of  time  was  obtained  for  each  test 
condition. 

To  evaluate  how  well  similarity  parameters  characterize  the  amount  of  suspended  sediment, 
experiments  were  repeated  over  a  range  of  rotor  speeds,  representing  different  wake  velocities, 
blade  passing  frequencies,  and  vortex  strengths.  Figure  2.33  shows  the  average  number  of  uplifted 
sediment  particles  for  each  recorded  rotor  revolution  as  a  function  of  rotor  rpm.  Each  data  point 
represents  the  average  number  of  particles  uplifted  over  one  rotor  revolution.  It  is  clear  that  the 
intensity  of  the  particle  field  during  each  revolution  varies  significantly,  even  for  a  given  rotor  rpm. 
For  the  steel  and  glass  sediment  particles  (Fig.  33(a)  and  Fig.  33(b)),  the  particle  field  intensity 
increased  with  rotor  speed. 

The  Buckingham-II  method  of  dimensional  analysis  was  employed  to  extract  the  relevant  scal¬ 
ing  parameters,  i.e.,  non-dimensional  groupings  of  parameters.  As  expected,  the  classical  similar¬ 
ity  parameters  arose  through  this  analysis,  but  so  did  new  parameters.  Based  on  previous  findings, 
the  vortices  are  known  to  affect  sediment  transport,  an  effect  not  captured  in  the  classical  analysis. 
To  characterize  this  critical  component  of  the  flow,  tip  vortex  characteristics  were  included  in  the 
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Figure  2.33:  Relative  brownout  intensity  at  different  rotor  speeds. 


dimensional  analysis  along  with  the  more  traditional  characterizations  of  the  rotor,  particles  and 
the  working  fluid.  Specifically,  the  rotor  was  characterized  by  its  radius,  chord,  tip  speed,  and  its 
rotational  frequency;  the  sediment  was  characterized  by  particle  diameter,  density,  terminal  veloc¬ 
ity  and  threshold  friction  velocity;  the  carrier  fluid  was  characterized  by  density,  a  characteristic 
velocity,  kinematic  viscosity,  and  boundary  layer  thickness  at  the  ground.  To  introduce  the  un¬ 
steady  aspects  of  the  flow  to  the  analysis,  the  tip  vortices  were  described  by  their  strength,  core 
size,  maximum  swirl  velocity,  and  area. 

The  resulting  parameters  were  paired  down  to  fewer  groupings.  Because  the  intention  was  that 
these  groupings  would  also  characterize  the  unsteady  features  of  the  flow,  it  was  required  that  the 
parameters  have  at  least  one  of  the  four  parameters  that  can  be  used  to  describe  the  vortices.  This 
requirement  eliminates  5  of  the  16  groupings.  Likewise,  two  groupings  were  eliminated  because 
they  did  not  include  sediment  properties.  Several  of  the  remaining  groupings  are  related  to  each 
other  by  one  of  the  six  classical  similarity  parameters;  for  example,  is  related  to  DlrUh  by 
Finally,  four  groupings  were  selected  and  one  additional  grouping  was  formed  by  the  product  of 
two  of  these.  The  final  candidates  were: 
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4.  —jrp f-,  called  the  threshold- swirl  velocity  ratio. 

Vq 

5.  — called  the  terminal/threshold- swirl  velocity  ratio. 

UpU* 

To  study  the  relationships  between  these  similarity  parameters  and  the  intensity  of  the  particle 
mobilization  and  uplift,  data  were  collected  across  a  range  of  rotor  operating  conditions  and  parti¬ 
cle  species.  The  terminal  velocity  and  threshold  friction  velocity  of  the  particles  was  varied  by  the 
selection  of  appropriate  particles.  The  vortex  strength  and  maximum  swirl  velocity  were  varied  by 
changing  rotor  speed. 

Figure  2.34  shows  results  for  the  quantity  of  particle  uplift  in  terms  of  each  of  the  five  candidate 
similarity  parameters.  Cases  with  multiple  data  points  are  experiments  that  were  repeated.  In 
Fig.  34(a),  the  quantity  of  particles  uplifted  is  plotted  with  respect  to  the  stationary  inertia  ratio 
and  shows  a  general  increasing  trend  for  each  sediment  type.  However,  this  similarity  parameter 
does  not  seem  to  properly  account  for  the  characteristics  of  the  sediment.  Figure  34(b)  shows  the 
same  results  in  terms  of  the  mobile  inertia  ratio.  Again,  each  sediment  type  shows  the  same  trend 
as  seen  in  Fig.  2.33  but  they  are  not  consistent.  The  same  outcome  was  obtained  for  remaining 
three  parameters.  The  threshold-swirl  velocity  ratio  results  shown  in  Fig.  34(d)  is  notable  as  it  is 
the  only  parameter  of  these  five  for  which  the  steel  and  glass  particle  results  overlap. 
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averaged  over  one  data  set  (1 000  images)  averaged  over  one  data  set  (1 000  images) 


(a)  Stationary  inertia  ratio  (b)  Mobile  inertia  ratio 


(c)  Terminal- swirl  velocity  ratio  (d)  Threshold-swirl  velocity  ratio 


(e)  Terminal-threshold-swirl  velocity  ratio 


Figure  2.34:  Quantity  of  particles  uplifted  in  terms  of  the  similarity  parameters. 
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Turbulent  Near-Bed  Flows  Under  Impinging  Vorticity 

Investigator(s):  J.  Gordon  Leishman 

Institution/Department:  University  of  Maryland/Department  of  Aerospace  Engineering 
Graduate  Student(s):  Jurgen  Rauleder 
Contact  email(s):  leishman@umd.edu 

Background  and  Objectives 

There  is  still  a  relatively  poor  understanding  as  to  how  discrete  turbulence  events  and  secondary 
flow  structures  affect  the  mobilization,  transport,  and  entrainment  of  particles  in  vortically  dom¬ 
inated  flows,  which  are  characteristic  of  those  found  in  brownout  conditions.  Also,  turbulence 
effects  have  not  yet  been  modeled  satisfactorily  for  such  flows.  Specifically,  for  the  vortically 
dominated  flows  found  in  the  wake  below  a  rotor,  time-resolved  measurements  of  the  detailed  tur¬ 
bulence  characteristics  are  not  available,  even  for  the  single-phase  flow  environment  (i.e.,  without 
dust  particles).  However,  detailed  measurements  for  such  a  flow  problem  are  needed  to  validate 
computational  models,  both  for  the  single-phase  and  also  for  the  dual-phase  dusty  flow.  Further¬ 
more,  the  sediment  entrainment  models  used  in  brownout  simulations  are  semi-empirical  in  nature 
and  have  been  derived  for  more  canonical  flows  other  than  rotor  flows  (i.e.,  steady,  uniform  bound¬ 
ary  layer  flows  for  the  most  part).  Therefore,  these  entrainment  models  do  not  explicitly  account 
for  unsteadiness  and  turbulence  in  the  flow,  and  they  have  not  been  validated  for  the  rotor-induced 
vortical  dual-phase  flow  involved  in  the  brownout  problem. 

The  pickup  and  suspension  of  particles  from  the  sediment  bed  is  at  least  partly  caused  by  suf¬ 
ficiently  high  surface  shear  stresses.  However,  the  actual  process  is  also  affected  by  fluctuations 
in  the  flow  velocities  and  the  turbulent  Reynolds  stresses,  discrete  turbulence  events,  and  the  un¬ 
steady  upward  flow  velocities  induced  by  the  vortices.  Most  prior  particle-laden  two-phase  flow 
studies  have  been  performed  on  channel  flows  with  relatively  “clean”  boundary  layers.  However, 
turbulence  effects  (including  the  effects  of  discrete  turbulent  motions)  and  the  potential  effects  of 
secondary  vortical  structures  on  particle  uplift  and  entrainment  have  not  been  studied  in  detail  in 
regard  to  the  underlying  flow-physical  mechanisms  in  non-uniform  flows. 

In  previous  work  under  this  task,  it  was  shown  that  it  was  the  rotor-generated  vortices  that 
were  primarily  responsible  for  sediment  uplift  from  an  underlying  sediment  bed.  It  was  observed 
that  particle  pickup  and  suspension  occurred  mostly  in  phase  with  the  periodic  passage  of  the  tip 
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vortices  over  the  bed,  but  the  role  of  turbulence  in  this  process  was  less  clear.  In  this  latest  research 
effort,  a  more  detailed  investigation  on  the  instantaneous  dual-phase  flow  field  at  the  ground  below 
the  rotor  has  been  carried  out.  Time-resolved  dual-phase  measurements  of  the  detailed  vortex  flow- 
particle  interactions  were  analyzed  to  further  the  understanding  of  the  fluid  dynamic  mechanisms 
involved  in  such  rotor-induced  particle  flows. 

Particular  attention  has  been  given  to  the  effects  of  organized  turbulent  motions  (or  turbu¬ 
lence  events)  and  secondary  vortical  structures  at  the  ground,  the  goal  being  to  better  understand 
the  fluid-particle  interactions  and  the  mechanisms  behind  particle  uplift.  The  quadrant  analysis 
method  was  also  used  to  expose  any  correlations  between  the  fluid  motion  and  the  particle  motion 
near  the  sediment  bed.  The  longer-term  objectives  were  to  help  with  the  development  of  more 
appropriate  dual-phase  flow  models  that  are  based  on  justifiable  physical  assumptions  that  could 
be  used  to  improve  brownout  simulations. 

Technical  Approach 

A  two-bladed  rotor  system  with  a  rotor  radius  of  0.085  m  (0.279  ft)  was  used  for  the  experiments. 
The  rotor  blades  used  thin  airfoils  (with  camber  and  a  sharpened  leading  edge)  of  rectangular 
planform  with  a  chord  of  0.018  m  (0.059  ft).  For  the  present  measurements,  the  rotational  plane 
of  the  rotor  was  located  at  a  height  of  one  rotor  radius  above  a  ground  plane  on  which  a  mobile 
sediment  bed  was  placed;  see  Fig.  2.35.  The  rotor  was  operated  at  a  rotational  frequency  of  60  Hz 
(3,600  rpm),  which  corresponded  to  a  blade  tip  speed  of  32.04  m  s-1  (105.12  ft  s-1)  and  a  chord 
Reynolds  number,  Rec,  at  the  blade  tip  of  approximately  35,000. 

The  blade  pitch  was  set  to  12°  to  produce  a  blade  loading  coefficient,  Cr/o,  of  0.156,  which 
yielded  a  theoretical  hover  induced  velocity,  v/j,  of  3.3  m  s-1  (10.8  ft  s-1).  From  the  measured 
velocity  field,  the  circulation  of  the  tip  vortices  was  determined  to  be  Tv  =  0.0172  m2s_1  (±  8.6  x 
10-4  m2s_1)  giving  a  vortex  Reynolds  number  of  Rev  =  1,100.  These  rotor  operating  conditions 
were  chosen  to  generate  a  representative  rotor  wake  with  enough  intensity  to  mobilize  and  uplift 
particles  of  sufficient  size,  to  allow  for  good  discrimination  of  the  flow  phases  in  the  resulting  dual¬ 
phase  flow,  as  well  as  to  avoid  significant  depletion  of  the  sediment  bed  during  the  measurement 
period. 

A  mobile  sediment  bed  consisting  of  well-characterized  glass  microspheres  with  a  diameter  of 
45-63  /am  was  used.  The  particles  were  loosely  deposited  on  the  ground  plane  to  an  approximate 
thickness  of  1  cm  and  the  surface  made  smooth  with  a  scraper.  The  Stokes  number  was  estimated 
to  be  60  and  the  particle  Reynolds  number  was  less  than  30.  The  simultaneous  time-resolved  flow 
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Figure  2.35:  Schematic  of  the  experimental  setup  used  to  conduct  the  dual-phase  PIV/PTV 
flow  measurements,  also  showing  the  coordinate  system  used  in  the  analysis. 

measurements  were  performed  using  two-dimensional,  two-component  PIV  for  the  carrier  phase, 
while  PTV  was  used  to  measure  the  motion  of  the  sediment  particles. 

A  high-repetition  rate  Nd:YLF  laser  and  a  4  Mp  CMOS  camera  (2,560  x  1,600  pixels  at  725 
frames  per  second)  were  used  to  interrogate  the  flow,  the  high  frame  rate  capability  allowing  for 
a  contiguous  time-history  of  the  flow  to  be  recorded.  This  time-history  enabled  detailed  assess¬ 
ments  as  to  how  the  fluid  flow,  its  primary  and  secondary  vortical  flow  structures,  and  turbulence 
events  near  the  sediment  bed  affected  the  particle  motion,  most  importantly,  how  they  affected  the 
processes  of  particle  uplift  and  entrainment.  The  selective  results  shown  in  this  report  are  pla¬ 
nar,  simultaneous  PIV/PTV  measurements  performed  in  the  regions  of  interest  (ROI)  as  shown 
in  Figs.  2.35  and  2.36.  An  exemplar  instantaneous  dual-phase  flow  measurement  below  the  rotor 
inside  the  larger  ROI  1  is  shown  in  Fig.  2.37. 

Results 

Although  the  younger  tip  vortices  at  radial  positions  of  r/R  ~  1.2-1. 3  convected  at  the  same  height 
above  the  sediment  bed  and  they  were  of  the  same  strength  or  stronger  than  the  older  vortices  (be¬ 
cause  of  vortex  distortion  and  diffusion  becoming  more  prevalent  further  downstream),  they  were 
mostly  unable  to  lift  up  particles  from  the  bed  at  these  upstream  locations;  see  Fig.  2.38.  There¬ 
fore,  in  this  reporting  period,  attention  of  the  present  task  was  given  to  secondary  flow  structures 
and  organized  turbulent  motions  as  they  may  contribute  to  the  complex  flow  environment  near  the 
bed  and  to  the  processes  of  sediment  mobilization  and  uplift.  To  visualize  the  turbulent  motions, 
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Figure  2.36:  Schematic  showing  the  regions  of  interest  (ROI)  that  were  investigated. 
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Figure  2.37:  Instantaneous  dual-phase  flow  realization  below  the  rotor.  Background  contour 
based  on  instantaneous  vorticity. 

it  is  necessary  to  observe  the  flow  in  a  moving  reference  frame.  Therefore,  the  mean  convection 
velocity  at  each  spatial  point  was  subtracted  a  posteriori  from  the  instantaneous  fluid  velocity  com¬ 
ponents  measured  by  the  PIV.  This  process  leaves  the  turbulent  fluctuations  about  the  mean  at  each 
location  in  the  flow  (i.e.,  a  Reynolds  type  of  decomposition),  as  shown  in  Fig.  2.38. 

For  the  instance  shown  in  Fig.  2.38,  the  tip  vortices  are  still  coherent  and  they  constitute  by 
far  the  strongest  vortical  structure  in  this  flow,  their  local  induced  velocities  being  much  greater 
than  the  local  mean  flow  velocities;  see  also  Fig.  2.37.  Therefore,  the  tip  vortices  are  apparent  in 
both  vector  maps,  i.e.,  both  when  the  flow  vectors  are  shown  on  basis  of  the  instantaneous  velocity 
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(a)  Instantaneous  velocity  vectors 


(b)  Instantaneous  fluctuation  velocity  vectors  (Reynolds  decomposed) 


Figure  2.38:  Velocity  vector  maps  superimposed  by  detected  particle  locations  for  the  same 
instance  in  ROI  2.  Only  every  other  measured  vector  shown  and  particle  sizes  exaggerated 
for  clarity. 


components  (Fig.  2.38(a))  and  when  shown  on  basis  of  the  Reynolds-decomposed  velocity  fluc¬ 
tuations,  u'  and  v';  see  Fig.  2.38(b).  However,  most  of  the  other  vortical  structures  and  organized 
turbulent  motions  in  the  flow  are  omitted  in  the  former  representation  of  the  flow,  as  shown  in 
Fig.  2.38(a).  The  Reynolds-decomposed  fluctuation  (or  perturbation)  velocity  gives  further  insight 
into  the  turbulent  flow  near  the  ground,  and  a  secondary  vortical  structure  becomes  visible  that  is 
counter-rotating  with  respect  to  the  tip  vortices  (vortical  center  at  r/R  =  1.4  and  z/R  =  0.12  in  this 
case). 

Notice  that  the  vortex  that  is  located  further  away  from  the  rotor  (viscous  core  at  r/R  =1.6 
in  this  case)  lifts  up  a  relatively  large  quantity  of  particles,  whereas  the  younger  tip  vortex  further 
upstream  does  not.  In  general,  the  younger  vortices  carried  more  concentrated  vorticity  than  the 
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vortices  further  downstream  because  they  were  more  coherent  and  less  diffused  than  the  vortices 
that  had  aged  more  in  the  flow  and  that  have  undergone  significant  interactions  with  the  sediment 
particles.  This  observation  indicated  that  the  tip  vortex  alone  was  not  always  able  to  lift  up  particles 
from  the  sediment  bed,  but  that  there  were  other  physical  mechanisms  that  can  also  contribute  to 
the  process.  From  Fig.  2.38(b),  notice  the  turbulent  motion  just  downstream  of  the  older  vortex 
(r/R  =  1.6-1. 8)  in  the  region  where  most  particles  are  uplifted  by  the  vortex.  This  decelerated  flow 
region  (in  the  streamwise  direction)  appears  to  convect  particles  onto  the  powerful  tip  vortex  that 
is  then  able  to  pick  up  the  particles  and  elevate  them  to  greater  heights  above  the  bed  and  even  into 
the  main  flow  closer  to  the  rotor. 

Figure  2.39  shows  consecutive  instantaneous,  Reynolds-decomposed  flow  realizations  near  the 
ground  (within  ROI  2)  for  every  30°  of  wake  age.  The  coherent  blade  tip  vortices  exposed  a  dis¬ 
tinct  quadrupole  structure  when  the  mV  correlation  was  plotted;  see  Fig.  2.39.  Most  tip  vortices 
trailed  a  decelerated  flow  region  (so-called  “inward  interaction”  motion,  for  r/R  =  1.6-1. 8  in  this 
case),  that  seemed  to  convect  particles  back  towards  the  succeeding  vortex,  where  they  were  subse¬ 
quently  picked  up  and  elevated  to  greater  heights  off  the  bed.  This  interesting  observation  rises  the 
question  if  this  turbulent  fluid  motion  was  a  precursor  to  particle  uplift,  and  if  it  was  a  necessary 
but  insufficient  condition. 

Figures  2.38(b)  and  2.39  suggested  that  positive  values  of  v'  (i.e.,  excess  velocities  in  the  up¬ 
ward  direction  away  from  the  bed)  prevailed  in  the  regions  where  most  of  the  sediment  particles 
were  picked  up  and  suspended  by  the  tip  vortices  (statistical  distributions  of  these  fluctuation  ve¬ 
locities  are  shown  later).  Furthermore,  an  ejection  motion  was  observed  downstream  of  most  tip 
vortices  and  it  contributed  to  the  process  of  particle  uplift.  This  ejection  motion  was  mostly  caused 
by  a  clockwise-rotating  secondary  vortical  structure  that  was  downstream  of  the  tip  vortex. 

The  quadrant  analysis  method  can  further  characterize  the  turbulence  based  on  organized  tur¬ 
bulence  structures  that  are  present  in  turbulent  flows  as  it  shows  the  correlative  structure  of  the 
fluctuation  (or  perturbation)  velocities,  u'  and  v' .  The  four  types  of  turbulence  events  (sweeps, 
ejections,  outward  and  inward  interactions)  are  classified  by  the  signs  of  the  fluctuation  velocities. 
The  measurements  are  then  plotted  on  a  quadrant  map  that  shows  the  probability  with  which  a 
certain  turbulence  event  occurs  at  one  spatial  location  in  the  flow.  This  quadrant  map  essentially 
shows  the  main  directions  of  the  turbulent  momentum  transfer  relative  to  the  bulk  movement  of 
the  fluid. 

Relatively  high  wall-normal  velocity  fluctuations,  v',  were  observed  throughout  all  downstream 
measurement  locations  at  a  height  of  z/R  =  0.1  above  the  sediment  bed;  see  Fig  2.40.  These 
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(c)  C  ~  720°  (d)  C  ~  750° 


Figure  2.39:  Consecutive  carrier  velocity  vector  fields  (Reynolds  decomposed)  superimposed 
by  the  particle  distribution  on  a  background  contour  showing  the  instantaneous  u'v'  correla¬ 
tion  in  ROI 2.  Only  every  other  measured  velocity  vector  is  shown. 

wall-normal  velocity  fluctuations  were  largely  produced  by  the  blade  tip  vortices.  Following  the 
development  of  the  flow  in  downstream  direction  for  this  constant  height  above  the  bed  shows 
a  shift  from  the  prevailing  streamwise  decelerated  fluid  ( r/R  <  1.4)  to  ejection  motions  (r/R  = 
1.5-1.65),  and  to  more  and  more  turbulence  events  with  greater  positive  wall-normal  fluctuations 
(v'  >  0;  upward,  away  from  the  ground)  at  r/R  =  1.65-1.8.  These  positive  (upwards)  excess 
velocities  are  important  for  particle  uplift  and  suspension  and  they  were  measured  to  be  more 
prominent  further  downstream  from  the  rotor  ( r/R  =  1.65-1.8),  which  correlated  well  with  the 
region  where  much  sediment  was  uplifted  to  sufficient  heights  above  the  ground  to  be  entrained 
into  the  main  vortex  flow;  see  Fig.  2.39. 

The  outward  interaction  events  gained  increasingly  more  weight  (i.e.,  their  occurrence  was 
more  probable)  downstream  of  r/R  =  1.6;  see  Fig.  2.40.  The  instantaneous  flow  vectors  based  on 
the  (Reynolds-decomposed)  velocity  fluctuations  shown  in  Figs.  2.39(c)  and  (d)  further  confirmed 
that  these  turbulence  events  were  induced  by  the  blade  tip  vortices  when  they  passed  by  and  lifted 


105 


Probability 


0  0.005 


0.01 


0.02  0.025  0.03 


(a)  r/R  =1.4 


(b)  r/R  =1.5 


(c)  r/R  =  1.6 


u\  u’\  uVvh 

(d)  r/R  =  1.65  (e)  r/R  =1.7  (f)  r/R  =  1.8 


Figure  2.40:  Joint  frequency  distributions  of  turbulence  events  for  the  carrier  at  z/R  =  0.1 
above  the  sediment  bed. 

up  particles  from  the  sediment  bed  in  this  region  ( r/R  ~  1.6-1. 8). 

Figure  40(d)  shows  the  statistically  prevailing  ejections  and  inward  interactions  around  r/R  ~ 
1.65  that  were  shown  to  facilitate  initial  particle  uplift;  see  e.g.,  Figs.  2.38  and  2.39.  These  decel¬ 
erated  ( u '  <  0)  turbulence  events  needed  to  be  followed  by  a  tip  vortex  convecting  downstream. 
In  most  cases,  only  with  this  interplay  between  the  decelerated  turbulent  motion  downstream  of 
the  vortex  and  the  positive  excess  velocities  induced  by  this  tip  vortex  itself,  may  the  particles 
then  be  uplifted  and  entrained  into  the  main  rotor  flow  above  the  sediment  bed,  as  also  previously 
discussed  in  Fig.  2.39. 

The  local  probability  density  distributions  for  the  turbulence  fluctuation  velocities  will  be  use¬ 
ful  for  modeling  and  simulation,  in  that  existing  numerical  models  can  be  validated.  Furthermore, 
future  computational  models  for  the  flow  near  the  ground  below  rotors  may  be  developed  by  in- 


106 


eluding  the  statistical  distributions  of  the  turbulence  fluctuations  given  by  this  task.  For  example,  in 
computational  models  that  are  unable  to  predict  turbulent  fluctuations  (i.e.,  models  that  are  merely 
descriptions  of  the  mean  flow),  these  perturbations  could  be  prescribed  by  statistical  probability 
distributions  to  better  simulate  the  flow  near  the  sediment  bed,  which  then  will  have  effects  on 
particle  mobilization  and  uplift. 
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Part  3:  Brownout  Synthesis,  Mitigation  and  Validation 


•  Task  3.1  -  Pilot-in-the-Loop  Brownout  Simulation  and  Mitigation  Strategies 

•  Task  3.2  -  Development  of  Efficient  Airborne  Sediment  Tracking  Algorithms 

•  Task  3.3  -  Understanding  Brownout  and  Developing  Mitigation  and  Control  Strategies 

•  Task  3.4  -  ABATE  Simulation  Framework  and  Validation 

•  Task  3.5  -  Computational  Considerations  in  the  Prediction  of  Brownout  Dust  Clouds 
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Task  3.1 


Pilot-in-the-Loop  Brownout  Simulation  and  Mitigation  Strategies 

Investigator(s):  Monica  Syal,  Chengjian  He,  Jin- Young  Hong,  and  Jinggen  Zhao 
Institution/Department:  Advanced  Rotorcraft  Technology,  Inc. 

Graduate  Student(s):  None 

Contact  email (s):  syal@flightlab.com,  he@flightlab.com 

Background  and  Technical  Challenges 

Modeling  brownout  and  the  associated  optical  degradations  in  a  flight  simulation  is  very  challeng¬ 
ing.  This  goal  requires  modeling  several  complex  physical  phenomena  in  real-time,  which  include: 
(a)  a  high-resolution  flow  field  produced  by  a  rotor  operating  in  ground  effect,  (b)  the  mobilization 
and  convection  of  a  very  large  number  of  dust  particles  over  relatively  long  time  scales,  (c)  the  op¬ 
tical  attenuation  caused  by  the  presence  of  dust  particles,  and  (d)  the  rendering  of  the  dust  clouds 
in  a  flight  simulator  environment.  Conducting  brownout  simulations  in  real-time  in  flight  simu¬ 
lators  requires  that  the  calculations  be  performed  in  less  than  20  ms  at  each  simulation  time  step. 
This  computational  goal  can  only  be  achieved  by  using  advanced  parallel  computing  techniques 
developed  from  scientific  computing,  and  by  tapping  the  computing  power  available  from  General 
Purpose  Graphics  Processing  Units  (GPGPUs). 

GPU-accelerated  brownout  simulation  methodology  was  developed  by  coupling  ART’s  flight 
dynamic  simulation  method  FLIGHTLAB,  Viscous  Vortex  Particle  Method  (VPM)  and  Particle 
Transfer  Method  (PTM,  developed  in  MURI  Task  3.2).  Recently,  the  focus  was  on  the  development 
of  physics-based  methodologies  to  “visualize”  brownout  dust  clouds  in  a  flight  simulator,  which 
includes  modeling  of  optical  obscuration  and  scattering  effects  of  dust  clouds  to  visualize  them  on 
graphic  processors.  Existing  image  generation  software  with  brownout  visualization  capabilities 
have  limited  applicability  because  they  simply  produce  brownout  dust  clouds  as  “rough  and  ready” 
special  effects  instead  of  representing  the  true  physical  aspects  of  the  problem,  i.e.,  the  correct 
optical  attenuation  of  light  based  on  spatial  frequency  and  the  correct  motion  of  the  dust.  To 
produce  realistic  and  physics-based  brownout  conditions  from  an  image  generator,  the  brownout 
simulation  methodology  needs  to  be  extended  to  render  the  optical  obscuration  and  scattering 
effects  of  the  dust  using  graphical  libraries  like  Open  Graphics  Library  (OpenGL). 
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Figure  3.1:  A  flow  chart  showing  integration  of  different  methodologies  to  develop  a 
brownout  simulation  tool  for  flight  simulators. 

Development  of  Methodologies 

Several  physics-based  methodologies  were  developed  to  calculate:  (a)  the  attenuation  of  light 
caused  by  brownout  dust  clouds,  (b)  the  rendering  of  dust  clouds  on  computer  graphics,  and  (c) 
the  integration  of  all  the  methodologies  to  simulate  and  render  dust  clouds  with  flight  simulator 
developed  at  Advanced  Rotorcraft  Technology  (ART)  Inc.  A  flowchart  showing  an  integration  of 
all  of  the  relevant  methodologies  is  shown  in  Fig.  3.1.  A  brief  discussion  of  the  new  methodologies 
developed  and  their  integration  is  given  below. 
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Figure  3.2:  A  flow  chart  showing  extraction  of  POD  mode  shape  <I>  and  coefficient  Q  matrices 
from  a  VPM  calculated  velocity  field  V. 

Rotor  Flow  Field  Using  Proper  Orthogonal  Decomposition  (POD) 

One  of  the  major  challenges  being  encountered  in  the  development  of  brownout  simulation  method¬ 
ology  for  flight  simulator  applications  is  the  real-time  modeling  of  rotor  wake  solution  using  VPM. 
The  cost  of  conducting  VPM  simulations  is  0(N {?),  where  Nv  is  the  number  of  vortex  particles  and 
for  a  good  fidelity  solution  Nv  ~  0(1O5).  While  efforts  are  underway  to  make  the  VPM  computa¬ 
tionally  faster  to  achieve  real-time  speedups,  an  interim  approach  has  been  followed  in  the  present 
work  to  achieve  real-time  performance.  In  this  approach,  a  pre-computed  VPM  solution  of  rotor 
wake  is  used  to  calculate  induced  velocities  at  dust  particles  in  real-time.  A  Proper  Orthogonal 
Decomposition  (POD)  [1]  based  algorithm  was  used  to  capture  the  dominant  features  of  the  VPM 
calculated  flow  field  using  only  a  very  small  number  of  “modes.”  This  technique  decomposes  the 
initial  velocity  data  set  into  separate  space  and  time  dependent  modes,  also  called  as  mode  shape 
functions  <E>(x)  and  coefficients  Q{t),  respectively;  see  Fig.  3.2.  In  the  present  work,  the  method¬ 
ology  to  reconstruct  induced  velocity  \p  on  the  dust  particle  positions  Xp  at  each  simulation  time 
step  using  POD  modes  was  implemented  in  parallel  on  GPUs  using  Compute  Unified  Device  Ar¬ 
chitecture  (CUDA)  to  achieve  real-time  performance;  see  Fig.  3.1. 

Development  of  Dust  Cloud  Illumination  (DCI)  Methodology 

The  appearance  of  a  dust  cloud  from  pilot’s  perspective  depends  upon  the  interaction  of  light  with 
the  dust  particles,  which  determines  the  “transparency  or  transmittance”  and  overall  “brightness” 
of  the  dust  cloud.  The  transmittance  T  of  a  medium  is  the  ratio  of  light  leaving  a  particular  point 
along  the  light  direction  and  reaching  another  point.  Brightness  (or  radiance)  I  of  the  dust  cloud 
determines  its  color  contrast  with  respect  to  the  surroundings.  The  computational  modeling  of 
these  effects  is  a  complex  problem.  In  the  present  work,  the  radiative  characteristics  of  light  and 
its  interactions  with  particles  were  modeled  using  the  Radiative  Transfer  Theory  (RTT)  [2], 

When  a  ray  of  light  interacts  with  matter  (e.g.,  particles),  it’s  electromagnetic  power  changes 
because  of:  (a)  absorption  (that  converts  the  light  energy  into  other  forms  of  energy,  e.g.,  heat 
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Figure  3.3:  A  schematic  showing  the  processes  of  scattering  and  transmission  of  the  incident 
light  along  direction  er  as  it  passes  through  a  cloud.  Source:  Harris  [2]. 

energy),  and  (b)  scattering  (that  changes  the  direction  of  light  after  it  interacts  with  the  matter). 
The  changes  in  the  power  of  electromagnetic  radiation  as  it  travels  through  the  dust  cloud  can  be 
calculated  using  RTT,  which  solves  for  the  changes  in  radiance  of  light  as  it  travels  through  the 
dust  cloud  because  of  the  effects  of  scattering  and  absorption. 

The  calculation  of  the  RTT  equation  can  be  explained  using  Fig.  3.3,  which  shows  a  ray  of  light 
incident  from  the  sun  along  a  direction  er  on  the  dust  cloud  as  perceived  by  an  observer.  The  RTT 
equation  evaluates  how  radiance  I  of  the  incident  light  changes  along  a  direction  er  (see  Fig.  3.3) 
as  it  travels  through  the  dust  cloud. 

RTT  accounts  for  two  types  of  effects:  (a)  single-scattering  effects  that  account  for  the  amount 
of  incident  radiance  that  gets  transmitted  through  the  dust  cloud  as  a  result  of  scattering  and  absorp¬ 
tion,  and  (b)  multiple-scattering  effects  that  account  for  the  changes  in  radiance  along  the  incident 
light  direction  because  of  scattering  from  multiple  particles  present  in  the  dust  cloud.  The  mod¬ 
eling  of  the  multiple-scattering  term  is  complex  because  it  involves  an  integration  of  the  scattered 
light  from  all  the  scattered  directions  from  all  the  particles  to  the  direction  of  the  incident  light  and 
is,  therefore,  usually  ignored  in  the  existing  models  to  predict  optical  attenuation  by  brownout  dust 
clouds  [3,  4],  However,  by  not  modeling  this  term  it  can  lead  to  artificially  darker  dust  clouds.  In 
the  present  work,  the  multiple-scattering  term  is  modeled  using  the  “multiple  forward  scattering” 
approximation,  in  which  multiple-scattering  of  light  is  only  considered  in  the  direction  of  the  eye. 
This  assumption  is  based  on  work  by  Harris  &  Lastra  [5], 

The  RTT  equation  has  been  solved  in  the  spherical  reference  frame.  The  particle  positions  are 
first  transformed  from  the  inertial  Cartesian  system  to  the  pilot-fixed  spherical  coordinate  system 
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using  the  methodology  developed  earlier  in  MURI  Task  3.1.  The  pilot’s  FOV  is  divided  into 
different  “bins,”  each  bin  being  characterized  by  its  solid  angle  and  its  radial  distance  from  the 
pilot’s  eye.  The  RTT  equation  is  solved  along  each  solid  angle  using  a  line  integral  method  (see 
details  in  Ref.  [2])  starting  from  outside  of  the  dust  cloud  and  moving  along  each  bin  in  the  solid 
angle  toward  the  pilot’s  eye.  The  solution  of  the  RTT  equation  is  transmittance  T  and  contrast  I  at 
each  bin.  These  parameters  T  and  I  are  used  to  determine  opacity  and  color  blending  for  rendering 
of  dust  clouds  in  flight  simulators.  In  the  present  work,  the  solution  of  the  RTT  equation  and  the 
calculation  of  T  and  I  have  also  been  fully  implemented  on  the  GPUs. 

Development  of  Dust  Cloud  Rendering  Methodology 

The  transmittance  and  brightness  of  dust  clouds  calculated  using  DCI  are  used  to  “render”  dust 
clouds  in  flight  simulators.  Rendering  means  creating  2D  images  from  3D  objects  using  computer 
graphics.  This  process  involves  a  camera  view  of  the  scene  based  on  a  light  source.  The  scene 
comprises  of  objects  (e.g.,  particles,  helicopter)  and  terrain,  which  are  converted  into  3D  triangles. 
Based  on  the  light  source,  these  triangles  defining  the  scene  are  then  converted  into  an  image  on 
the  computer  screen. 

The  dust  clouds  are  rendered  on  computer  graphics  using  an  OpenGL  based  open-source  soft¬ 
ware  visualization  toolkit  OpenSceneGraph  (OSG).  This  toolkit  is  written  in  standard  C++  and 
OpenGL,  and  can  run  on  most  of  the  operating  systems.  In  the  present  work,  a  particle  render¬ 
ing  software  was  developed  based  on  in-built  OSG  data  structures  by  defining  the  physical  and 
graphical  properties  of  particles.  The  physical  parameters  include  number  Np,  positions  Xp  and 
velocities  \p  of  particles  provided  by  PTM  and  VPM.  The  graphical  or  optical  parameters  include 
size,  contrast  I  and  transparency  T  (calculated  using  DCI). 

Integration  with  Flight  Simulator 

The  OSG-based  rendering  methodology  developed  in  this  work  was  integrated  as  a  plugin  with 
an  open-source,  cross  platform  image  generator  called  Multi  Purpose  Viewer  (MPV).  Figure  3.1 
shows  that  all  of  the  methodologies  used  in  this  work  have  been  coupled  into  two  separate  pro¬ 
cesses,  which  are  launched  in  parallel  in  the  FLIGHTLAB  simulation  environment  using  two  CPU 
threads  in  the  same  computer.  In  the  first  process,  FLIGHTLAB,  VPM  solution  (using  POD  ma¬ 
trices),  PTM  and  DCI  methodologies  are  linked  together  to  calculate  particle  velocity  \p,  position 
Xp,  transparency  T,  and  contrast  I.  All  of  these  methodologies  are  integrated  by  passing  data  point- 
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ers  on  the  GPU  memory,  i.e.,  all  methods  are  integrated  within  the  GPU  memory  using  CUDA.  At 
the  beginning  of  the  simulation,  all  of  these  methodologies  are  initialized  and  the  data  is  transferred 
to  the  GPU  memory.  There  after,  all  of  these  methodologies  interchange  data  with  each  other  with 
in  the  GPU  memory.  This  approach  was  followed  to  reduce  expensive  CPU-GPU  memory  trans¬ 
fers.  The  outcomes  of  the  first  process,  i.e.,  \p,  Xp,  T,  and  I  are  integrated  to  a  second  process, 
which  is  running  MPV  with  the  brownout  rendering  software  developed  in  this  work  as  a  plugin. 
This  integration  is  performed  using  a  shared  memory  block  on  the  computer’s  main  memory. 

Results 

1.  Rendering  of  Dust  Clouds: 

Dust  clouds  are  rendered  (in  OSG)  with  graphical  properties,  i.e.,  transparency  and  contrast, 
calculated  using  the  Dust  Cloud  Illumination  (DCI)  methodology.  The  results  showing  ren¬ 
dered  dust  clouds  were  obtained  for  a  representative  single-rotor  helicopter  in  an  approach 
maneuver  developed  from  Moen  et  al.  [6].  The  approach  simulation  was  conducted  for 
about  1 1  seconds.  A  total  of  8.0 x  105  monodisperse  particles  with  a  diameter  of  20  /./in  were 
used  in  this  study. 

The  pilot’s  FOV  was  divided  into  300,000  bins.  The  transparency  and  radiance  of  light  that 
reaches  pilot’s  eye  depend  on  particle  concentration  r\  in  each  bin,  which  is  the  number  of 
particles  per  unit  volume  of  each  bin.  Typically,  the  value  of  T|  varies  between  0(1O10)  and 
0(  101 )  for  particles  of  sizes  dp  varying  between  0(1)  /am  and  0(1000)  /mi,  respectively  [7, 
8].  Representing  such  high  concentration  of  particles  in  dust  cloud  simulations  is  impractical 
because  of  the  limitations  of  computer  memory  and  the  required  simulation  time.  Therefore, 
in  the  present  work  particle  concentrations  were  scaled  to  represent  high  concentration  for 
the  calculations  of  transparency  and  brightness  (or  radiance). 

Figures  4(a),  4(b)  and  4(c)  show  the  instantaneous  realizations  of  three  rendered  dust  clouds 
(using  OSG)  in  the  pilot’s  FOV  with  particle  concentrations  O(106),  O(108)  and  O(109), 
respectively.  These  figures  show  that  as  T|  increases,  a  dust  cloud  becomes  less  transparent 
because  of  the  presence  of  more  particles  in  the  optical  path.  The  dust  clouds  are  rendered 
in  the  present  work  by  rendering  individual  particles.  There  is  a  paradox  in  this  approach 
because  the  transparency  and  brightness  are  scaled  to  achieve  r|  O(108)  to  0(1O9).  How¬ 
ever,  rendering  such  a  large  number  of  particles  in  real-time  becomes  impractical.  One  way 
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(c)  r|  =  0(1  ()9j  (d)  r|  =  0(109),  bigger  sized  particles 


Figure  3.4:  Instantaneous  realizations  of  dust  clouds  at  t  =  1 1  s  obtained  using  particle  con¬ 
centrations:  (a)  T|  =  0(1 06) ,  (b)  0(1 08 ),  (c)  0(1 09 ) ,  and  (d)  0(1 09 )  with  bigger  size  particles 
to  represent  high  particle  concentrations. 

to  achieve  the  effect  of  such  large  number  of  particles  is  by  representing  each  particle  by  a 
bigger  texture  that  can  cover  more  number  of  pixels,  see  Fig.  4(d). 

Figures  5(a)  and  5(b)  show  instantaneous  realizations  of  dust  clouds  in  the  OpenSceneGraph 
(OSG)  by  using  only  the  single-scattering  term,  and  including  the  multiple-scattering  term 
to  calculate  brightness  of  dust  clouds,  respectively.  It  is  apparent  from  the  two  figures  that 
the  dust  cloud  visualized  using  the  multiple-scattering  model  is  significantly  brighter  when 
compared  to  the  single-scattering  models  because  multiple-scattering  produces  higher  con¬ 
trast  in  the  forward  direction  because  of  scattering  from  all  of  the  particles  in  the  dust  cloud. 
Therefore,  it  is  important  to  model  multiple-scattering  in  calculating  light  attenuation. 


2.  Integration  of  All  Methodologies  with  ART’s  Flight  Simulator: 

NVIDIA’s  Kepler  K10  GPU  was  used  in  the  present  work.  Figure  3.6  shows  a  variation  of  the 
taken  taken  by  the  K10  GPU  to  calculate:  (a)  the  induced  velocity  using  the  pre-calculated 
POD  matrices,  (b)  the  particle  positions  using  PTM,  and  (c)  the  transparency  and  contrast 
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(a)  Radiance  calculated  using  single-scattering  (b)  Radiance  calculated  using  multiple-scattering 
model  model 


Figure  3.5:  Instantaneous  realizations  of  a  dust  cloud  with  rendering  calculated  using:  (a) 
only  single-scattering  model,  and  (b)  both  single  and  multiple  scattering  models. 

using  DCI.  The  number  of  bins  used  in  DCI  to  bin  pilot’s  FOV  are  300,000.  Notice  that  for 
NP  <  106,  all  of  these  simulations  collectively  require  less  than  10  ms  to  be  calculated  on 
one  GPU.  For  the  current  goal  of  Np  «  105,  these  simulations  require  <  3  ms,  i.e.,  they  can 
be  calculated  in  real-time. 

Figure  3.7  shows  set-up  of  the  flight  simulator  in  ART  with  the  brownout  simulation  plugin. 
Instantaneous  realizations  of  brownout  dust  clouds  obtained  during  an  approach  maneuver  in 
the  flight  simulator  are  shown  in  Fig.  3.8.  These  figures  show  that  the  brownout  simulation 
methodology  has  been  successfully  integrated  with  the  flight  simulator. 
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Figure  3.6:  Variation  of  average  time  taken  by  GPU  to  execute  different  methodologies,  i.e., 
reconstruction  of  velocity  field  using  POD  modes,  PTM  and  DCI. 


Figure  3.7:  Set-up  of  ART’s  flight  simulator  integrated  with  brownout  plugin. 


117 


(e)  (f) 


Figure  3.8:  Dust  cloud  images  of  a  brownout  simulation  captured  from  ART’s  flight  simula¬ 
tor. 
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Task  3.2 


Development  of  Numerically  Efficient  Airborne  Sediment  Tracking  Algorithms 

Investigator(s):  J.  Gordon  Leishman,  Nail  Gumerov 

Institution/Department:  University  of  Maryland/Department  of  Aerospace  Engineering 
Graduate  Student(s):  Monica  Syal,  Bharath  Govindarajan 
Contact  email:  leishman@umd.edu 

Background  and  Technical  Challenges 

A  significant  issue  in  realistic  dust  cloud  simulations  is  the  very  large  number  of  particles  that 
are  needed  to  obtain  dust  clouds  of  acceptable  fidelity.  Computing  the  motion  of  each  and  every 
individual  sediment  particle  in  a  dust  cloud  is  a  computationally  very  expensive  problem,  and  so 
novel  solution  approaches  must  be  used.  Additional  technical  challenges  include  the  modeling  of 
the  detailed  processes  of  particle  mobilization  and  lift-off  from  a  sediment  bed,  the  convection  of 
the  particles  in  a  complex  vortical  flow,  the  potential  for  particle-particle  interactions  if  the  dust 
clouds  become  sufficiently  dense,  and  the  effects  of  previously  suspended  particles  colliding  with 
the  sediment  bed,  i.e.,  more  particles  that  are  liberated  by  bombardment  mechanisms. 

Technical  Approach 

This  task  involves  the  development  of  computationally  efficient  algorithms  that  can  be  applied  to 
the  simulation  of  dilute  gas-particle  suspensions  at  low  Reynolds  numbers  of  the  relative  particle 
motion.  Computationally  expeditious  algorithms  have  been  developed  to  solve  for  the  particle 
dynamics  in  a  Lagrangian  framework,  including  both  the  uplift  from  the  sediment  bed  and  the  sub¬ 
sequent  convection  through  the  flow.  Surface  shear  stresses  and  unsteady  pressure  forces  mostly 
determine  the  threshold  for  the  mobilization  of  particles  from  the  bed. 

After  the  particles  are  lofted  into  the  flow  from  the  bed,  their  trajectories  are  computed  us¬ 
ing  particle  equations  of  motion  (EOMs)  in  the  Lagrangian  frame  of  reference.  These  three- 
dimensional  equations  are  a  coupled  set  of  ordinary  differential  equations  (ODEs),  involving  iner¬ 
tial  terms  (gravity,  centrifugal,  Coriolis),  and  aerodynamic  terms  (drag  and  lift).  Efficient,  stable 
and  time-accurate  multi-step  predictor-corrector  based  numerical  integration  schemes  have  been 
explored  to  solve  the  particle  EOMs.  Of  course,  the  fidelity  of  the  approach  depends  upon  the 
number  of  dust  particles  that  can  be  simulated  in  reasonable  time  frames.  Because  it  may  re¬ 
quire  0(1O13)-0(1O16)  individual  particles  to  produce  a  realistic  brownout  cloud,  the  resulting 
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approaches  have  been  parallelized  to  run  on  multiple  processors  as  well  as  Graphical  Processing 
Units  (GPUs)  (see  Task  3.5),  and  also  takes  advantage  of  particle  clustering  schemes. 


Particle  Mobility  Modeling 

Stationary  particles  on  a  bed  below  the  rotor  can  experience  several  types  of  forces  such  as 
shear,  pressure,  inter-particle,  and  gravitational.  Shear  forces  on  the  stationary  particles  arise  be¬ 
cause  of  the  velocity  gradients  produced  by  the  boundary  layer  on  the  sediment  bed.  Significant 
unsteady  pressures  may  also  be  produced  in  the  vorticity-laden  flow  field  below  the  rotor.  Besides 
these  forces,  a  particle  on  the  bed  also  experiences  gravitational  forces  and  inter-particle  forces. 
A  new  model  has  been  developed  in  this  task  to  represent  the  process  of  particle  mobilization  by 
considering  the  effects  of  both  shear  and  pressure  forces.  The  expression  of  the  threshold  friction 
velocity,  uH ,  which  is  the  minimum  friction  velocity  required  to  mobilize  particles,  is 


(3.1) 


where  AP  is  the  unsteady  pressure  difference  acting  on  the  particle,  pp  is  particle  density,  dp  is 
particle  diameter,  p  is  flow  density,  g  is  the  acceleration  under  gravity,  A  is  a  threshold  friction 
velocity  parameter  that  depends  upon  particle  Reynolds  number,  and  y  is  a  parameter  that  depends 
upon  the  inter-particle  forces  acting  on  the  particle. 

Sediment  Tracking 

Once  the  particles  are  uplifted  into  the  main  flow,  the  primary  forces  acting  on  the  particles 
are  aerodynamic  drag  and  gravitational  forces.  The  resulting  particle  motion  is  governed  by  the 
Basset-Boussinesq-Oseen  (BBO)  equations,  which  under  the  assumptions  of  dilute  gas-particle 
suspensions  of  low  Reynolds  number  Stokes’  flow,  can  be  decoupled  in  the  three  spatial  dimen¬ 
sions.  The  resulting  equations  becomes  two  sets  of  three  first-order  ODEs,  i.e., 


(3.2) 


(3.3) 


where  i  =  x,  y,  z,  and  f  is  the  forcing  function  containing  the  contributions  from  the  flow  velocity, 
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V,  and  the  gravitational  force,  and  can  be  written  as 

fx  =  Vx 
fy  =  Vy 

fz  =  Vz-gxPs  (3.4) 

where  xPs  is  the  particle  response  time  given  by  xPs  =  ppdp/lS/j. 

Several  time-marching  integration  schemes  have  been  developed  in  this  task  to  solve  the  par¬ 
ticle  equations  of  motion  (EOMs),  as  given  by  Eqs.  3.2  and  3.3.  These  schemes  were  developed 
by  considering  their  stability,  accuracy  and  convergence  characteristics  for  the  practical  time  step 
sizes  that  would  be  needed  in  a  brownout  simulation.  Most  of  the  explicit  numerical  integration 
schemes  are  unstable  for  solving  such  stiff  equations  unless  a  very  small  time  step  size  is  used. 
The  stability  of  implicit  schemes,  such  as  Euler  implicit,  two-point  and  three-point  backward  dif¬ 
ference  formulae,  were  examined.  Several  predictor-corrector  based  time-marching  algorithms 
using  the  BDF2  and  BDF3  schemes  were  also  developed,  with  better  results,  but  they  incur  a 
higher  cost  because  two  induced  velocity  field  calculations  are  needed.  The  accuracy  of  these 
schemes  have  been  analyzed  for  several  prototypical  problems,  and  the  the  schemes  based  on  the 
BDF2  and  BDF3  methods  were  second-  and  third-order  accurate,  respectively.  A  time-marching 
method  based  on  the  BDF2  scheme  was  ultimately  used  for  the  dust  cloud  simulations  because  its 
accuracy  is  consistent  with  the  accuracy  of  the  free- vortex  wake  solution. 

Bombardment  Ejections 

A  suspended  particle  can  also  impact  the  ground  and  so  lose  part  of  its  kinetic  energy.  De¬ 
pending  upon  its  momentum  before  the  impact,  as  well  as  the  inelasticity  and  frictional  forces  pro¬ 
duced  by  the  sediment  bed,  the  particle  will  rebound.  Some  of  this  energy  goes  into  ejecting  more 
particles  from  the  bed,  which  is  called  bombardment  ejections;  see  Fig.  3.9.  The  bombardment 
mechanism  was  modeled  by  assuming  that  the  particle  impacting  the  bed  excavates  void-shaped 
craters,  the  volume  of  particles  emitted  being  equal  to  the  volume  of  the  crater  formed  by  the  im¬ 
pacting  particle.  The  initial  locations  of  the  ejected  particles  is  based  on  a  probabilistic  approach 
using  energy  conservation.  Both  calculations  and  experiments  have  shown  that  bombardment  ejec¬ 
tion  is  one  of  the  primary  modes  of  sediment  uplift  during  the  development  of  brownout  conditions. 

Use  of  GPU  Programming 

One  of  the  primary  challenges  in  the  modeling  of  brownout  dust  clouds  is  the  tracking  of  the 
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Impacting  particle  Ejected  particles 


Figure  3.9:  Schematic  of  a  sediment  particle  impacting  the  ground,  forming  a  crater,  and 
ejecting  new  particles. 

typically  very  large  number  of  particles  over  the  time  scales  characteristic  of  the  flight  maneuvers, 
i.e.,  when  a  rotorcraft  is  performing  a  landing  or  takeoff.  To  contain  this  cost,  both  the  aerodynamic 
and  dust  cloud  simulations  were  structured  so  that  they  could  be  executed  on  multiple  processors. 
The  particle  simulation  algorithm  was  also  implemented  on  graphic  processor  units  (GPUs)  by 
dividing  and  assigning  particle  data  to  each  GPU  thread  and  spawning  a  large  number  of  similar 
threads.  Further  details  of  this  approach  are  given  in  Task  3.5. 

Clustering  Methods 

Particle  clustering  methods  have  been  examined  in  detail  for  their  potential  application  to  the 
problem  of  predicting  rotorcraft  brownout.  Computationally  efficient  algorithms  were  developed 
that  can  be  applied  to  the  simulation  of  dilute  gas-particle  suspensions.  In  particular,  the  Gaussian 
distribution,  k-means  and  Osiptsov’s  clustering  methods  were  examined  for  a  prototypical  flow 
field  that  mimics  the  highly  unsteady,  two-phase  vortical  particle  flow  found  during  brownout 
conditions. 

The  Gaussian  method  was  generally  the  most  successful  approach  to  clustering.  In  this  method, 
the  individual  particles  in  the  brownout  simulation  are  treated  as  clusters,  each  representing  a  spe¬ 
cific  number  of  particles.  The  clusters  are  convected  at  each  time  using  the  EOMs.  The  positions 
of  the  individual  particles  associated  with  a  cluster  are  calculated  by  using  a  multi- variate  Gaussian 
distribution.  The  position  of  the  particles  within  a  specific  cluster  are  recomputed  at  every  time, 
and  therefore,  can  be  implemented  as  a  post-processing  step  in  the  simulation  framework. 

Figure  3.10  shows  an  example  of  the  application  of  the  Gaussian  clustering  technique  to  a  three 
dimensional  brownout  cloud.  A  comparison  of  the  solution  was  obtained  with  106  cluster  centers 
(denoted  by  the  black  dots  in  Fig.  3.10(a))  and  using  the  distributed  particle  clusters  (Fig.  3.10(b)) 
with  a  total  of  100  particles  per  cluster,  giving  a  total  of  108  particles.  Notice  that  the  clustered 
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(D)  Zoomed  in  view  to  show, 
details 


Details  of 
computed  dust 
cloud  are  always 
maintained 


(A)  Actual  simulated  solution  / 


(B)  Clustered  (brown)  and  actual 
solution  (black) 


(C)  Particle  densities  from  the 
clustered  solution 


Figure  3.10:  A  solution  obtained  from  brownout  dust  field  computations  when  using  the 
Gaussian  clustering  method. 

solution,  as  shown  by  Fig.  3.10(b),  identifies  the  actually  computed  cluster  centers  as  black  dots, 
and  the  particles  obtained  through  the  Gaussian  clustering  technique  as  brown  dots.  The  results 
clearly  show  that  the  clustered  solution  brings  out  the  finer,  structured  details  of  the  dust  cloud. 
Figure  3.10(c)  shows  the  density  of  the  dust  cloud  that  is  derived  from  the  clustered  solution,  and 
Fig.  3.10(d)  shows  some  details  of  the  dust  cloud  near  the  ground,  the  concentrations  of  which  can 
be  correlated  to  the  positions  of  the  tip  vortices  in  the  flow. 

Two  other  clustering  methods  have  also  been  explored.  The  k-means  clustering  method  is  based 
on  the  principle  that  certain  sets  of  individual  particles  can  be  decomposed  into  smaller  groups  of 
clusters,  and  that  the  resulting  EOMs  are  solved  only  for  the  clusters.  There  are  two  aspects  to 
this  method:  1.  Selecting  the  candidate  particle  groups  to  form  a  cluster;  2.  The  effects  produced 
by  clustering,  declustering,  and  reclustering  of  the  particle  groups.  The  third  method,  Osiptsov’s 
method,  has  its  roots  in  full  Lagrangian  theory  and  involves  integrating  the  equations  for  particle 
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density  along  individual  pathlines.  The  particle  densities  are  obtained  from  the  Lagrangian  form 
of  the  conservation  of  mass  by  computing  the  change  in  volume  of  an  element  of  “particle  fluid” 
along  its  pathline. 
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Understanding  Brownout  and  Developing  Mitigation  and  Control  Strategies 

Investigator(s):  Roberto  Celi,  J.  Gordon  Leishman,  and  J.  Baeder 
Institution/Department:  University  of  Maryland,  Department  of  Aerospace  Engineering 
Graduate  Student(s):  John  Tritschler,  Jillian  Alfred 

Contact  email(s):  celi@umd.edu,  leishman@umd.edu,  baeder@umd.edu,  chopra@umd.edu 

Background  and  Technical  Challenges 

The  development  of  a  brownout  dust  cloud  is  a  complex  phenomenon  that  is  affected  by  many 
parameters  that  define  the  characteristics  of  the  rotorcraft  as  well  as  by  the  flight  profile  and  landing 
trajectory  that  is  followed  by  the  pilot.  The  sediment  bed  itself  also  introduces  complexities  into  the 
brownout  problem  in  that  it  has  properties  and  characteristics  that  vary  from  location  to  location. 

The  mitigation  and  control  of  brownout  is  a  correspondingly  complex  problem:  (i)  “Mitiga¬ 
tion”  must  be  defined  through  some  quantitative  metric  or  metrics,  (ii)  the  design  methodology 
must  capture  the  fundamental  physics  of  all  potential  mitigation  solutions,  (iii)  numerical  opti¬ 
mization  is  an  ideal  tool  to  quantify  and  achieve  mitigation,  but  it  is  difficult  to  apply  because  of 
the  large  computational  requirements  associated  with  brownout  modeling  as  well  as  the  potential 
complexities  of  the  design  space  (e.g.,  its  nonconvexity,  time  dependency,  disjointness).  Finally, 
brownout  mitigation  cannot  be  the  sole  design  objective,  and  additional  requirements  (e.g.,  vehicle 
performance,  vibration  levels,  handling  qualities,  etc.)  need  to  be  considered. 

Technical  Approach 

In  this  task,  the  problem  of  brownout  has  been  better  understood  by  running  simulations  that 
expose  the  underlying  mechanisms  of  dust  cloud  formation  around  a  helicopter  rotor  that  is  in 
dynamic  motion,  i.e.,  one  that  is  moving  relative  to  the  ground  plane.  The  brownout  analysis  used 
is  that  described  under  Task  3.2.  Dust  particles  that  are  originally  stationary  on  the  ground  are 
mobilized  by  the  action  of  the  rotor  wake,  mainly  by  the  shear  stresses  induced  on  the  sediment 
bed  by  the  tip  vortices,  which  also  cause  the  particles  to  be  convected  upwards  and  then  into 
suspension.  As  the  concentration  of  particles  increases  in  front  of  the  rotor,  the  pilot’s  visibility 
of  the  ground  and  the  landing  zone  will  diminish,  which  in  practice  leads  to  a  loss  of  situational 
awareness. 

This  research  has  examined  the  problem  of  brownout  mitigation  from  several  perspectives:  1. 
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Mitigation  through  flight  path  management;  2.  Mitigation  by  means  of  changes  to  the  rotor  design; 
3.  Mitigation  by  means  of  direct  changes  to  the  characteristics  of  the  blade  tip  vortices,  i.e.,  by 
diffusing  the  vortices.  In  each  case,  the  mitigation  activities  have  been  formulated  as  a  formal 
numerical  optimization  problem.  The  objective  function  B(X)  that  has  been  used  in  the  present 
research  is  based  on  the  space  and  time  distributions  of  the  particle  concentrations,  weighing  more 
heavily  areas  of  good  visibility  that  allow  the  pilot  to  close  attitude  and  position  loops,  i.e.,  the 
areas  immediately  in  front  and  to  the  side  of  the  pilots.  While  this  metric  is  only  a  surrogate  for 
the  actual  optical  attenuation  of  the  light  as  it  passes  through  the  brownout  cloud,  this  objective 
function  has  proved  successful  in  the  present  work  in  the  sense  that  its  minimization  has  led  to  the 
prediction  of  dust  clouds  around  the  rotor  that  are  less  concentrated  in  the  pilot’s  field  of  view,  i.e., 
that  there  are  fewer  particles  uplifted  there. 

Results 

To  mitigate  the  severity  of  brownout,  previous  work  performed  under  this  task  has  examined  how 
several  aspects  of  the  rotor  design  (including  the  effects  of  rotor  radius,  blade  chord,  blade  twist 
rate,  and  the  number  of  blades)  can  affect  brownout  severity.  Results  were  obtained  for  a  4- 
bladed  rotor  that  is  representative  of  a  medium  size  utility  helicopter  of  a  prescribed  weight  that 
is  executing  and  approach  to  landing  over  a  sediment  bed.  It  has  been  shown  that  for  a  given 
number  of  blades  the  reducing  the  disk  loading  resulted  in  reduced  brownout  severity.  At  the 
same  time,  the  range  of  different  results  that  are  obtained  from  the  same  (or  similar)  disk  loading 
values  indicates  that  brownout  severity  cannot  be  attributed  to  the  effects  of  disk  loading  alone. 
The  research  has  also  shown  that  the  severity  of  brownout  is  closely  associated  with  the  blade 
tip  vortex  strengths.  For  example,  the  combined  effects  of  simultaneously  increasing  number  of 
blades  and  reducing  the  individual  tip  vortex  strengths  leads  to  greater  brownout  severity,  which 
arises  from  the  complicated  effects  that  these  parameters  have  on  the  interactions  of  the  tip  vortices 
with  the  ground  and  with  each  other. 

The  recent  activities  performed  under  this  task  have  been  related  to  brownout  mitigation  by 
means  of  flight  path  management.  The  previous  brownout  mitigation  study,  which  consisted  of 
the  determination  of  optimum  landing  trajectories,  was  repeated  with  a  more  sophisticated  simu¬ 
lation  model.  Whereas  the  previous  model  was  limited  to  an  isolated  rotor  with  rigid  blades  and 
prescribed  blade  control  angles,  the  new  model  consisted  of  a  full  flight  dynamic  simulation,  in¬ 
cluding  a  detailed  representation  of  rotor  blade  dynamics.  The  brownout  analysis  (see  Task  3.2) 
was  coupled  with  a  comprehensive  flight  dynamic  model,  which  included  first  principles  descrip- 
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tions  of  aircraft  and  rotor  dynamics.  In  last  year’s  mitigation  study,  the  thrust  on  the  isolated 
rotor  and  the  corresponding  controls  along  the  trajectory,  were  determined  through  a  procedure 
resembling  trim  at  every  time  step.  The  values  of  the  rotor  thrust  were  plausible  but  not  exact.  In 
the  latest  model,  both  the  controls  and  the  rotor  forces  and  moments  were  those  required  to  fol¬ 
low  the  desired  trajectories.  Aircraft  body  motions  and  blade  flapping  responses  were  also  those 
appropriate  to  the  flight  trajectories. 

The  basic  trajectory  optimization  methodology  was  the  same  as  that  used  previously.  The 
objective  function  to  be  minimized  was  a  measure  of  the  particle  concentration  of  the  cloud  (i.e., 
particle  density  in  the  flow  in  terms  of  numbers  of  particles  per  unit  volume)  at  each  instant  of 
the  terminal  phase  of  landing.  Behavior  constraints  were  used  to  protect  the  rotor  from  entering 
the  vortex  ring  state  conditions  (an  undesirable  and  potentially  unsafe  flight  condition),  and  from 
entrance  into  dangerous  areas  of  the  height- velocity  diagram  where  a  safe  engine-out  autorotational 
landing  would  not  be  feasible.  The  solution  methodology  was  iterative.  At  each  optimization 
step,  first,  a  genetic  algorithm  was  used  to  find  all  the  minima  (the  optimization  problem  is  not 
convex)  of  a  response  surface-type  approximation  of  the  objective  function.  Next,  the  optimum 
was  analyzed  precisely,  i.e.,  a  more  thorough  brownout  simulation  was  carried  out.  Finally,  the 
precise  objective  function  was  used  to  update  the  response  surface  approximation. 

The  design  variables  used  in  the  optimization  were  the  same  as  that  used  in  the  previous  op¬ 
timization  study,  i.e.,  three  parameters  describing  the  flight  path  and  landing  trajectory.  It  was 
assumed  that  the  helicopter  was  equipped  with  a  flight  control  system  (FCS)  that,  for  each  design 
vector  selected  by  the  optimizer,  i.e.,  for  each  trajectory,  the  FCS  would  adjust  the  aircraft  controls 
so  as  to  follow  that  trajectory.  The  control  laws  were  obtained  using  a  simple  LQG/LTR  method¬ 
ology.  This  is  a  “pseudo-FCS”  actually  meant  to  represent  a  human  pilot,  albeit  an  especially  well 
trained  one,  and  one  that  is  not  affected  by  the  potential  loss  of  visual  cues.  This  optimization 
formulation  had  the  same  general  convergence  characteristics  and  computational  requirements  as 
that  solved  previously,  but  the  recent  results  are  considered  to  be  more  realistic. 

The  original  plans  for  current  study  called  for  the  formulation  of  the  optimization  as  an  inverse 
simulation  problem,  where  the  design  variables  would  have  been  the  time  histories  of  the  pilot 
controls.  This  situation  is  the  most  realistic  formulation  but  it  is  computationally  much  more 
intensive  (in  this  case  it  would  have  required  at  least  12-16  design  variables).  The  procedure 
described  above,  with  the  “pseudo-FCS”  was  intended  to  be  a  preliminary  step  to  the  inverse 
simulation,  to  provide  good  initial  guesses  of  the  controls  and  identify  the  regions  of  the  design 
space  with  local  minima.  Instead,  it  proved  to  be  a  very  effective  optimization  procedure  on  its 
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Figure  3.11:  Time  history  of  body  positions  of  the  actual  flight  trajectory  compared  to  the 
desired  approach  profile  for  X  =  [6°  1 10  kts  300  ft]r. 

own.  The  resulting  trajectories  for  one  of  the  cases  are  shown  in  Fig.  3.11.  The  simulated  pilot 
controls  were  also  very  realistic.  Some  spurious  high  frequency  content  was  present  in  some  phases 
of  the  maneuver,  but  these  had  negligible  effects  on  the  brownout  simulation. 

The  key  conclusions  of  the  recent  study  are: 

1.  It  is  possible  to  define  a  controller  that  will  generate  the  time  histories  of  the  controls  re¬ 
quired  to  track  a  specified  approach  profile,  in  place  of  a  more  rigorous  but  more  compu¬ 
tationally  intensive  inverse  simulation  procedure.  A  sufficiently  accurate  trajectory  resulted 
from  customizing  the  controller  for  the  flight  condition  at  each  time  step.  The  pilot  stick 
inputs  remained  within  a  range  that  was  realistic  relative  to  their  hard  limits. 

2.  The  results  of  the  study  with  a  full,  free-flight  model  of  the  helicopter,  largely  confirmed 
those  obtained  with  the  simpler,  isolated  rotor  only  model,  in  particular,  that  the  roll  up  of 
the  blade  tip  vortices  in  the  rotor  wake  and  the  formation  of  a  ground  vortex  ahead  of  the 
rotor  plays  a  key  role  in  the  formation  of  a  brownout  dust  cloud. 
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Figure  3.12:  Dust  cloud  geometries  top  view  for  (a)  the  baseline  approach  profile  and  (b)  the 
shallow  optimum  approach. 

3.  The  design  space  for  the  approach  profile  optimization  was  non-convex  and  contained  at 
least  two  optima:  a  local  optimum  corresponding  to  a  shallow  angle  approach  in  which  the 
pilot  “outruns”  the  cloud  that  is  forming  behind  the  helicopter  until  the  actual  landing,  and  a 
global  optimum  corresponding  to  a  steep  trajectory  in  which  the  pilot  effectively  keeps  the 
rotor  wake  further  above  the  ground  for  as  long  as  possible  until  landing.  These  results  are 
known  to  be  consistent  with  current  operational  practice,  but  there  are  also  other  possible 
variations.  Dust  cloud  geometries  for  the  baseline  and  the  optimum  steep  approach  profiles 
are  shown  in  Fig.  3.12. 

4.  While  the  optimum  approach  trajectories  were  found  to  be  very  similar  to  that  of  the  previous 
study,  which  was  based  on  a  simplified  rotor-only  model,  the  characteristics  of  the  resulting 
brownout  cloud  were  different.  The  addition  of  aircraft  dynamics  to  the  simulation  changed 
the  position  of  the  center  of  gravity  of  the  aircraft,  and  hence  changed  the  rotor  tip-path-plane 
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Figure  3.13:  Realizations  of  the  velocity  magnitudes  and  developing  cloud  in  a  longitudinal 
plane  through  the  flow  field  for  (a)-(e)  the  baseline,  Xi,  and  (f)-(j)  the  global  (steep)  optimum, 
X20,  approach  profiles.  The  motion  of  the  rotor  over  the  ground  is  from  right  to  left.  Particles 
are  enlarged  for  illustrative  purposes. 

angle  of  attack.  The  difference  in  blade  flapping  response  caused  the  blade  tip  vortices  to 
impinge  on  the  ground  further  inboard  of  the  blade  tips,  resulting  in  the  dust  particles  being 
mobilized  and  uplifted  closer  to  the  rotor.  The  forward  pitching  of  the  rotor  also  resulted  in 
greater  fuselage  pitch  up  angles  seen  during  the  maneuver.  Furthermore,  changing  the  center 
of  gravity  of  the  helicopter  changed  the  characteristics  of  the  resulting  brownout  cloud.  An 
example  of  the  velocity  magnitudes  and  developing  cloud  in  a  longitudinal  plane  through 
the  flow  field  for  the  global  optimum,  corresponding  to  a  steep  descent  approach,  is  shown 
in  Fig.  3.13. 
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ABATE  Simulation  Framework  and  Validation 

Investigator(s):  James  Baeder 

Institution/Department:  University  of  Maryland,  Department  of  Aerospace  Engineering 

Graduate  Student(s):  Sebastian  Thomas 

Contact  email (s):  baeder@umd.edu,  sthomas2@umd.edu 

Background  and  Technical  Challenges 

Rotorcraft  brownout  problem  is  a  complex,  multi-physics  phenomenon,  the  modeling  of  which, 
requires  the  capability  to  reliably  predict  the  rotor  loads  and  the  vortex  wake  generation  from 
the  rotor  in  ground  effect,  the  interaction  of  the  wake  with  the  fuselage,  and  the  dual-phase  flow 
field  near  the  ground  involving  the  rotor  wake  and  sediment  particles  entrained  into  the  flow  field. 
A  first-principles  based  continuum  dynamics  treatment  of  the  brownout  problem  is  impractical. 
However,  a  modular  approach  using  domain  decomposition  principles,  for  e.g.,  domain  decom¬ 
position  at  the  fluid- structure  interface,  is  a  more  practical  option.  Such  an  approach  allows  the 
individual  solvers  to  use  the  most  efficient,  domain- specific  solution  technique  to  solve  the  gov¬ 
erning  equations.  This  type  of  approach  also  allows  greater  flexibility  in  the  development,  testing 
and  validation  of  the  individual  solvers,  thereby  providing  greater  confidence  in  the  final  coupled 
simulation  framework  as  a  predictive  tool. 

The  design  of  a  robust,  scalable  and  numerically  accurate  modular  framework  poses  many 
technical  challenges.  In  addition  to  the  numerical  accuracy  of  the  participating  solvers,  care  must 
be  taken  to  ensure  that  the  coupling  processes  preserve  the  desired  order  of  accuracy  and  do  not 
introduce  additional  errors  during  information  transfer  across  domain  boundaries.  The  modular 
framework  must  have  a  well-defined  interface,  over  which  the  participating  solvers  can  exchange 
data  with  other  solvers  without  having  extensive  knowledge  of  the  implementation  details  of  the 
participating  solvers.  Such  capability  is  extremely  valuable  in  the  development  and  testing  of 
newer  solvers  by  permitting  easier  coupling  with  previously  validated  solvers.  Because  the  sim¬ 
ulation  of  the  brownout  problem  will  inevitably  involve  the  use  of  high-fidelity  CFD  solvers,  the 
capability  to  support  execution  in  a  parallel  environments  (using  domain  decomposition  strategies 
or  shared  memory  approaches)  in  a  transparent  manner  is  of  utmost  importance.  The  framework 
must  be  capable  of  managing  the  execution  of  participating  solvers  across  multiple  processors  or 
on  modern  platforms  like  Graphics  Processing  Units  (GPUs),  and  must  ensure  seamless  informa- 
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Figure  3.14:  Schematic  illustrating  the  interaction  of  the  various  MURI  tasks. 

tion  transfer  not  only  between  different  solvers,  but  also  across  different  processors  without  the 
intervention  from  participating  solvers  or  the  end  user. 


Technical  Approach 

The  primary  objective  of  this  task  is  to  synthesize  a  highly  efficient,  Advanced  Brownout  Aerome¬ 
chanic  Technology  Environment  (ABATE)  simulation.  The  approach  builds  upon  expertise  in  com¬ 
prehensive  rotorcraft  aeromechanics  analyses  while  incorporating  strategies  and  computational 
models  that  are  developed  and  used  in  the  other  MURI  tasks  -  RANS  solvers,  Lagrangian  vortex 
methods,  sediment  uplift  and  suspension  models,  as  well  as  optical  transmissibility  through  the 
resulting  dust  clouds  (Task  3.1).  Figure  3.14  summarizes  the  various  interactions  between  the 
different  MURI  tasks. 

The  primary  purpose  of  Task  3.4  is  the  integration  of  the  computational  models  that  are  devel¬ 
oped  and  validated  under  the  other  MURI  tasks.  For  the  simulation  of  each  individual  domain,  the 
framework  allows  the  user  to  choose  from  multiple  models  varying  in  sophistication  and  compu¬ 
tational  expense.  In  the  present  approach,  the  fluid  phase  is  modeled  using  a  hybrid  methodology 
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Figure  3.15:  Schematic  describing  the  structure  of  the  Hybrid  FVM-RANS  Solver.  (1)  For 
the  airloads  distribution,  a  linearized  aerodynamics  module  is  used.  (2)  Computations  in  the 
far-wake  region  between  the  rotor  tip  path  plane  are  performed  using  a  free- vortex  method. 
(3)  A  high-fidelity  RANS  solver  is  used  for  computations  near  the  ground  plane.  (4)  For 
dual-phase  flows,  the  RANS  solution  is  coupled  to  a  particle  transport  solver. 

combining  the  capabilities  of  a  RANS  solver  with  a  free-vortex  wake  method.  The  central  idea 
involves  the  usage  of  the  grid  based  solver  near  the  ground  plane  where  complex  phenomena  like 
vortex-ground  interactions  and  boundary  layer  formation/separation  occur.  The  precise  location 
and  clustering  of  this  ground  mesh  is  chosen  based  on  a  careful  analysis  of  the  full-RANS  results 
from  Task  1 .2.  For  regions  that  are  well  separated  from  the  ground,  a  filament  based  free-vortex 
method  (FVM),  similar  to  the  model  used  in  Task  3.2,  is  employed. 

Figure  3.15  illustrates  the  combination  of  RANS  and  free-vortex  wake  that  is  used  to  model 
the  fluid-phase.  At  the  RANS-FVM  interface,  a  non-reflecting  boundary  condition  is  employed  to 
admit  vortical  structures  into  the  Eulerian  mesh.  Figure  3.16  shows  a  representative  ground  mesh 
system  used  in  the  hybrid  methodology. 

For  the  dual-phase  calculations,  the  sediment  phase  is  currently  being  modeled  using  a  La- 
grangian  particle  tracking  algorithm  (see  Task  3.2),  which  is  one  way  coupled.  In  the  current 
work,  all  simulations  are  performed  in  a  time-accurate  fashion. 

This  hybrid  approach  leads  to  a  large  reduction  in  computational  expense  when  compared  to 
full  RANS  simulations  that  employ  Eulerian  grids  to  simulate  the  entire  domain.  While  this  strat¬ 
egy  is  superior  to  the  full-RANS  simulations  in  terms  of  numerical  efficiency,  for  peak  performance 
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(a)  Isometric  view  of  the  ground  mesh  used 
for  RANS  calculations.  The  dimensions  of 
this  mesh  are  60x180x120 
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(b)  Side  view  of  the  ground  mesh  showing 
refinement  to  resolve  both  tip  vortices  and  the 
boundary  layer 


Figure  3.16:  Mesh  system  used  in  the  hybrid  simulation  of  the  micro-rotor. 

it  is  still  required  to  run  the  hybrid  method  on  dozens  of  processors  using  a  domain  decomposition 
strategy.  For  a  given  problem  size,  domain  decomposition  is  not  scalable  to  an  arbitrarily  large 
number  of  processors  -  in  this  limit,  the  time  taken  for  information  transfer  between  sub-domains 
becomes  comparable  to  the  time  taken  for  computations.  To  circumvent  this  difficulty,  a  new 
approach  was  sought. 

In  Task  3.2,  one  of  the  strategies  employed  for  algorithmic  speedup  involved  the  use  of  GPU 
technology  to  accelerate  the  highly  data-parallel  Biot-Savart  computations  in  the  free-vortex  wake 
model.  The  resulting  solver  demonstrated  considerable  speedup  over  equivalent  serial  versions. 
Using  this  performance  gain  as  a  benchmark,  all  components  of  the  hybrid  method  were  rewritten 
and  optimized  in  CUDA-C.  Each  component  was  then  validated  using  canonical  test  cases  and 
benchmarked  on  various  gaming  and  general-purpose  GPU  platforms.  On  combining  these  dif¬ 
ferent  GPU-based  solvers,  the  resulting  desktop  simulations  using  the  hybrid  method  were  shown 
to  be  significantly  faster  than  previous  calculations  involving  domain  decomposition  strategies. 
Previously,  a  single  revolution  of  an  MAV-scale  rotor  required  close  to  6  hours  on  32  Intel  Xeon 
cores.  Now,  using  the  GPU-based  hybrid  method  on  a  gaming  platform  like  the  GTX  Titan,  a 
single  revolution  is  completed  in  just  1.25  hours. 
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(a)  Mesh  system  used  for  RANS  calculations  in  the  impinging  vortex  ring  simulations 
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(b)  Hybrid  methodology  prediction  of  vortex  ring  trajectory 

Figure  3.17:  CFD  validation  using  hybrid  methodology  of  experiments  performed  under 
Task  2.3. 

Results 

The  GPU-based  hybrid  methodology  was  first  validated  with  experiments  performed  under  task 
2.3.  The  experiments  consisted  of  an  impinging  wall  jet  formation  by  modulating  the  flow  at  the 
exit  plane  of  a  nozzle  to  produce  a  coherent  vortex  ring  that  proceeded  to  interact  with  the  ground 
plane.  Figure  3.15(a)  shows  the  CFD  mesh  system  that  consisted  of  a  cylindrical  background  mesh 
for  RANS  calculations.  The  free-vortex  wake  filaments  in  the  Lagrangian  wake  were  initialized  in 
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(a)  Formation  of  sediment  waves  seen  in  the  experiments 
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(b)  Formation  of  sediment  waves  seen  in  the  simulations 

Figure  3.18:  The  formation  of  sediment  waves  outboard  of  the  rotor.  Note:  The  particles  are 
not  drawn  to  scale. 

a  ring  topology  with  the  end  of  the  last  filament  with  the  beginning  of  the  first.  The  impinging  vor¬ 
tex  ring  was  modeled  using  an  analytical  expression.  Figure  3.15(b)  shows  the  tip  vortex  trajectory 
predicted  by  the  hybrid  methodology  by  superimposing  seven  representative  time  instances  in  the 
simulation.  The  formation  of  a  secondary  vortical  structure  was  observed  at  the  ground  plane  sim¬ 
ilar  to  a  phenomena  reported  in  experiments.  The  approximate  position,  size  and  strength  of  this 
secondary  vortex,  relative  to  the  primary  vortex  were  also  found  to  be  in  good  agreement  with  ex¬ 
perimental  observations.  The  computations  were  quantitatively  validated  by  comparing  predicted 
time-averaged,  velocity  profiles  at  the  ground  plane  with  experimental  measurements. 

The  hybrid  methodology  was  then  used  to  model  the  single-phase  and  dual-phase  how  envi¬ 
ronment  beneath  a  hovering  micro-scale  rotor.  The  simulations  were  performed  in  a  time-accurate 
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(a)  Isometric  view  of  the  dual-phase  flow  field  at  the  start  of  the  computations.  Both  the  free- 
vortex  wake  away  from  the  ground  and  iso-surfaces  of  vorticity  inside  the  RANS  mesh  are  plotted 
to  show  the  tip  vortex  trajectory.  Most  particles  are  lying  immobile  on  the  ground  plane. 


(b)  Isometric  view  of  the  dual-phase  flow  field  after  16.25  revolutions.  The  dispersed  phase  is 
highly  active  with  particles  being  entrained  around  close-proximity  vortices  to  form  well  defined 
three-dimensional  waves. 


Figure  3.19:  Predicted  3D  structure  of  the  dual-phase  flow  field  with  50  micron  particles 
beneath  a  hovering  MAV-rotor.  Note:  The  particles  are  not  drawn  to  scale. 

fashion.  The  single-phase  results  were  validated  with  experiments  performed  under  Task  1.1.  The 
work  under  Task  2.2  with  a  micro-scale  rotor  operating  in  ground  effect  over  a  sediment  bed  had 
identified  six  primary  transport  mechanisms  that  dominated  the  dispersed  phase  in  vortical  flows. 
These  experimental  results  were  used  to  validate  computational  approaches  that  simulate  the  two- 
phase  flow  field.  In  the  present  approach,  sediment  particles  were  mobilized  and  entrained  into 
the  flow  using  threshold  friction  velocity  condition  based  on  Bagnolds  model.  Unsteady  pres- 
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sure  effects  and  particle  collisions  were  neglected,  though  the  modular  structure  of  the  ABATE 
environment  means  that  other  sediment  mobility  models  can  be  used.  Figure  3.18(a)  compares  a 
snapshot  of  the  dual-phase  simulation  with  an  image  visualizing  the  particle  distribution  seen  in 
the  laboratory  experiments  (Fig.  3.18(b))  conducted  under  Task  2.2.  The  creation  of  sediment 
waves  was  observed  outboard  of  the  rotor  similar  to  those  seen  in  experiments,  if  multiple  layers 
of  particles  existed  in  the  sediment  bed.  Figure  3.19  shows  an  isometric  snapshot  of  a  dual-phase 
simulation  showing  the  3D  structure  of  a  typical  brownout  cloud.  Clearly  defined  and  sustained 
three-dimensional  sediment  waves  are  seen  to  form  in  simulations  involving  particles  of  50  mi¬ 
crons.  This  size  of  particles  were  also  observed  to  exhibit  the  mechanisms  of  creep  and  vortex 
trapping  with  the  latter  mechanism  being  the  dominant  mode  of  transport  at  this  scale.  Further¬ 
more,  the  regions  of  sediment  mobility  were  seen  to  coincide  with  the  radial  zones  where  the 
predicted  friction  velocity  was  observed  to  be  the  largest. 
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Task  3.5 


Computational  Considerations  in  the  Prediction  of  Brownout  Dust  Clouds 

Investigator(s):  Ramani  Duraiswami,  Nail  A.  Gumerov,  J.  Gordon  Leishman 
Institution/Department:  University  of  Maryland,  University  of  Maryland,  Department  of  Aerospace 
Engineering  &  Department  of  Computer  Science 
Graduate  Student(s):  Qi  Hu,  Monica  Syal 
Contact  email(s):  leishman@umd.edu 

Background  and  Technical  Challenges 

The  primary  goal  of  this  task  was  to  develop  and  test  high  performance  computing  (HPC)  methods 
capable  of  efficiently  handling  the  large-scale  numerical  problems  that  occur  in  the  computation 
of  rotorcraft  brownout.  There  are  several  approaches  that  can  be  used  to  compute  such  flows,  in¬ 
cluding  Eulerian  methods  based  on  finite  difference  discretization  and  Lagrangian  methods.  The 
latter  are  among  the  most  advanced  of  contemporary  computational  fluid  dynamics  methods.  The 
present  work  makes  use  of  Lagrangian  methods,  although  the  modeling  and  computational  ap¬ 
proaches  being  developed  in  this  research  task  are  not  restricted  to  only  these.  Lagrangian  methods 
are  mesh  free  (or  may  require  only  surface  meshing  as  opposed  to  volumetric  meshes),  but  require 
the  tracking  and  interaction  of  large  numbers  of  flow  singularities  and  particles,  or  particle  clusters. 

The  technical  challenges  faced  included:  1.  The  effective  utilization  of  advanced  hardware, 
including  single  and  multiple  graphics  processors  (GPUs)  as  well  as  clusters  consisting  of  many 
heterogeneous  nodes,  each  of  which  consists  of  multicore  CPUs  and  several  GPUs;  2.  The  de¬ 
velopment  and  use  of  advanced  scalable  algorithms  that  reduce  computational  cost  of  the  problem 
from  O  (N2  +  NM )  to  O  (N  +  M),  thus  providing  orders  of  magnitude  improvements  in  computa¬ 
tional  speed.  Such  algorithms  are  also  memory  compact,  which  is  very  important;  3.  Integrating 
the  cumulative  effects  of  hardware  and  algorithmic  accelerations,  i.e.,  the  implementation  of  the 
advanced  algorithms  on  advanced  heterogeneous  architectures,  which  requires  solution  of  many 
technical  tasks  such  as  heterogeneous  programming  and  the  modification  of  the  algorithms  for  the 
most  efficient  use  of  workstations  and  clusters,  etc.;  4.  Challenges  related  to  mathematical  mod¬ 
eling  to  reduce  the  computational  complexity  of  the  problems,  which  includes  brownout  particle 
clustering  algorithms,  developing  a  combination  of  methods  of  continuum  mechanics  of  multi¬ 
phase  flows  with  direct  simulations,  development  of  efficient  interpolation  techniques,  stabiliza¬ 
tion  of  numerical  schemes  (e.g.,  the  use  of  high-order  symplectic  time  marching  schemes,  error 
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control,  and  the  use  of  heterogeneous  precision  computations),  and  frameworks  for  comparisons 
with  other  numerical  techniques. 

Technical  Approach 

Modeling  the  brownout  problem  is  complex  and  requires  expensive  computations  of  two-phase, 
unsteady,  three-dimensional  vortical  flows.  The  memory  and  speed  limits  of  conventional  compu¬ 
tational  resources,  such  as  personal  computers  and  small  or  mid-size  CPU  clusters  present  a  barrier 
in  the  use  of  computations  to  conduct  parametric  studies,  solve  optimization  and  design  problems, 
and  investigate  the  brownout  phenomenon  based  on  first  principles  modeling.  The  research  con¬ 
ducted  in  this  task  has  focused  on  several  directions  to  help  solve  these  problems: 

1.  Hardware  acceleration.  The  use  of  graphics  processors  (GPUs)  for  acceleration  of  already 
developed  and  tested  methods  for  brownout  simulations; 

2.  Algorithmic  acceleration.  Performing  fundamental  research  on  the  development  of  a  scal¬ 
able  fast  multipole  method  (FMM)  on  heterogeneous  computing  architectures,  which  can 
substantially  accelerate  simulations  by  the  realization  of  scalable  algorithms; 

3.  Combined  hardware/algorithmic  acceleration.  The  design  of  efficient  data  structures  and 
homogeneous/heterogeneous  algorithms  for  single  PCs,  multicore  PCs,  multicore  PCs  with 
graphics  processors  and  heterogeneous  clusters,  which  have  computing  nodes  that  consist  of 
multicore  CPUs  and  several  GPUs. 

While  initial  work  in  items  2  and  3  was  focused  on  the  algorithms  themselves,  the  technology 
was  transferred  to  other  members  of  the  MURI  team.  Also,  with  a  view  to  extending  the  algorithms 
to  situations  involving  tail  rotors  and  fuselages,  research  was  conducted  on  fast-multipole  acceler¬ 
ated  vortex-lattice  methods  and  versions  for  these  new  algorithms  have  been  developed  and  tested. 
Special  attention  was  paid  to  the  development  of  a  user-friendly  Matlab  front  end  environment  for 
CPU/GPU  heterogeneous  programming.  Research  and  development  of  basic  algorithms  has  also 
continued,  so  that  they  can  be  accelerated  and  applied  to  large  scale  problems.  The  latter  includes 
extension  of  the  algorithms  to  problems  with  periodic  boundary  conditions. 

Results 

The  research  has  proceeded  toward:  1.  Validation  and  testing  of  the  methods,  and  2.  Finding 
solutions  to  bottleneck  problems  that  affect  the  performance  and  scalability  of  the  algorithms. 
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Figure  3.20:  Brute  force  acceleration  of  brownout  simulations  using  the  GPU. 


Graphics  Processors  and  Brute  Force  Acceleration 

This  activity  is  related  to  the  efficient  utilization  of  graphics  processors  to  accelerate  existing 
algorithms  and  methods.  The  parts  of  the  FORTRAN  algorithm  related  to  computation  of  the 
velocity  field  at  the  particle  location  and  translation  of  particles  were  mapped  on  the  GPU  using 
FLAGON  middleware  library  [1].  The  dust  clouds  during  a  landing  maneuver  were  computed  for  a 
representative  single  rotor  helicopter  with  a  4-bladed  rotor.  A  total  of  1  million  particles  were  used 
in  this  case  to  obtain  the  dust  cloud.  To  obtain  35  seconds  of  actual  simulation  time  (approximately 
150  rotor  revolutions)  of  the  dust  cloud  development  computed  in  a  CPU  only  environment  takes 
about  23  days  of  wall-clock  time.  If  implemented  on  the  GPUs,  the  same  solution  can  be  obtained 
in  23  hours,  and  without  compromising  on  accuracy  [2],  However,  the  actual  performance  gains 
obtained  from  the  GPU  depend  upon  the  problem  parameters  and  the  desired  numerical  precision. 

Figure  3.20  shows  the  computational  performance  gains  obtained  for  such  brownout  simula¬ 
tions,  implemented  on  both  single  precision  and  double-precision  GPU  architectures  for  different 
particle  counts.  The  performance  gain  is  defined  as  a  ratio  of  wall-clock  times  required  to  run  a 
case  on  the  GPU  architecture  and  the  full  CPU  simulation.  The  results  obtained  are  clearly  very 
good:  250  times  speed-up  for  single-precision  and  25  times  speed-up  for  double-precision.  No- 
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tice  that  the  latter  numbers  are  obtained  using  NVIDIA  Cl 060  cards,  for  which  there  is  almost  an 
order  of  magnitude  difference  in  single  and  double  precision  simulations.  Newer  cards,  such  as 
the  NVIDIA  C2050  have  reduced  this  gap  substantially  (only  two  times  difference  in  the  perfor¬ 
mance). 

Important  things  learned  during  this  task  were  related  to  the  stability  of  the  time  marching 
schemes  with  single  and  double  precision,  multi-language  programming,  and  interfacing  of  the 
conventional  FORTRAN  and  C/C++  codes  with  graphics  processors,  which  basic  programming 
was  performed  using  the  CUDA  environment  and  the  middleware  library  mentioned  previously. 

Algorithmic  Accelerations 

Methods  to  compute  high  Reynolds  number  vortical  flows  include  the  Eulerian,  Eulerian- 
Lagrangian,  and  fully  Lagrangian  methods.  Each  approach  has  its  relative  advantages.  The  present 
work  has  focused  on  Lagrangian  methods,  which  are  better  suited  to  handle  enormous  amount  of 
spatially  non-uniform  unstructured  data  because  they  avoid  large  volumetric  meshes.  The  compu¬ 
tation  of  single-phase  flows  can  be  undertaken  using  vortex  element  methods  (VEM),  while  the 
computation  of  the  two-phase  flow  can  be  realized  by  means  of  the  tracking  of  individual  repre¬ 
sentative  particles  in  a  given  carrier  fluid  flow  or  in  a  flow  affected  by  the  dust  particles  by  using 
one-  and  two-way  interaction  schemes,  respectively;  see  Task  3.2. 

The  major  challenge  in  the  efficient  use  of  Lagrangian  methods  in  the  simulation  of  brownout  is 
the  modeling  of  the  large  number  of  singularities  (e.g.,  markers  in  the  rotor  wake)  that  are  required 
to  represent  the  flow  with  the  necessary  level  of  detail  and  accuracy.  Each  Lagrangian  marker 
also  may  interact  with  another  Lagrangian  marker  (vortex  element  or  other  fluid  particle),  so  the 
number  of  such  interactions  at  each  time  step  becomes  N2,  where  N  is  the  number  of  markers. 
In  the  case  when  a  second  phase  (particle  or  dust)  is  present  and  is  represented  by  M  particles, 
there  are  an  additional  NM  interactions.  Therefore,  the  overall  computational  complexity  of  the 
problem  becomes  extremely  large,  i.e.,  O  (N2  +NMj.  In  terms  of  algorithm  developments  two 
major  improvements  can  be  envisioned.  First,  reduction  of  M  via  some  intelligent  scheme.  Second, 
the  use  of  scalable  algorithms  of  complexity  0(N+M )  or  with  some  additional  logarithmic  factors 
with  respect  to  N  and  M. 

The  first  problem  can  be  attacked  by  using  clustering  algorithms,  continuum  descriptions  of  the 
multiphase  flow,  and  fully  Lagrangian  approaches  for  the  particle  number  density.  Several  algo¬ 
rithms  of  these  type  have  been  tried  and  showed  some  improvement  [3],  but  technical  challenges 
still  remain,  especially  in  algorithms  related  to  particle  number  density  tracking  because  they  re- 
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quire  advanced  interpolation  techniques  for  unstructured  data  sets.  Several  algorithms  have  been 
developed  for  this  purpose,  however  they  still  require  more  study  for  high  performance  implemen¬ 
tations  and  so  far  comprehensive  tests  have  been  conducted  only  in  two  dimensions. 

One  of  the  major  algorithmic  developments  of  the  present  research  is  related  to  the  fast- 
multipole  accelerated  vortex-element  methods.  A  new  translation  theory  has  been  developed  for 
such  methods  based  on  the  Lamb-Helmholtz  decomposition  of  the  vortical  fields  [10].  Such  a  for¬ 
mulation  and  theory  enables  reduction  of  computations  of  velocity  and  velocity  gradient  fields  to 
the  evaluation  of  just  two  coupled  scalar  harmonic  potentials  (solutions  of  the  Laplace  equation). 
In  some  reported  high  performance  computing  of  vortical  flows  using  the  fast  multipole  methods, 
the  computation  of  the  vortex  stretching  terms  was  found  to  be  very  computationally  expensive. 
The  algorithms  developed  in  this  research  task  were  shown  to  reduce  the  computational  cost  by 
several  times. 

Another  algorithmic  improvement  was  the  efficient  computation  of  vortex  core  functions.  Sev¬ 
eral  versions  of  the  VEM  were  developed,  including  vortex  particle  method  and  vortex  filament 
method,  which  were  carefully  tested  on  a  number  of  well-known  problems,  such  as  the  interaction 
of  vortex  rings  and  filaments.  The  method  was  also  tested  for  some  experimental  cases  conducted 
in  Task  1.1.  One  such  comparison  is  presented  in  Fig.  3.21.  Here,  measurements  of  the  veloc¬ 
ity  field  generated  by  a  single  rotating  blade  that  were  obtained  using  particle  image  velocimetry 
(Task  1.1)  were  compared  to  results  using  the  vortex  method.  It  can  be  seen  that  the  measured  and 
simulated  fields  are  in  good  general  agreement. 

Fast  Multipole  Methods 

The  vortex  methods  are  accelerated  in  the  present  work  via  the  fast  multipole  methods  (FMM). 
Solution  of  practical  problems  using  a  Lagrangian  approach  with  scaling  0(N 2)  makes  it  impos¬ 
sible  to  simulate  problems  larger  than  N  ^  105.  The  FMM  is  the  break-through  algorithm,  which 
increases  this  limit  by  several  orders  of  magnitude.  One  of  the  major  contributions  of  this  project 
was  the  development  of  a  much  faster  distributed  parallel  version  of  the  FMM. 

In  the  beginning  of  this  task,  a  basic  commercial  FMM  code  was  available  for  the  scalar 
Laplace  equation  in  three  dimensions.  This  code  was  used  as  a  kernel  and  routines  were  built 
to  accommodate  the  new  translation  theory  mentioned  above.  In  addition,  new  methods  were 
developed  to  create  an  FMM  on  distributed  parallel  computing  architectures.  In  the  developed  al¬ 
gorithm  each  component  had  fine  grain  parallelization  for  multicore  machines,  while  the  work  was 
distributed  by  a  novel  strategy  across  many  machines.  The  whole  algorithm  was  revised,  profiled. 


145 


Simulation 


Experiment 

v, ,  m/s 


x,  mm 
v,,  m/s 


60  80  100  120  140  160  ISO 

x,  mm 

|v|,  m/s 


V,, ,  m/s 


60  80  100  120  140  160  180 

x,  mm 
v,,  m/s 


|v|,  m/s 


60  80  100  120  140  160  180 

x,  mm 


20 

80  1 

60 

15 

mm 

o 

10 

M 

r 

20 

o  1 

-10 


60  80  100  120  140  160  180 

x,  mm 


Figure  3.21:  Comparison  of  the  measured  and  simulated  velocity  field  components,  vx,  vz, 
and  |v||  |  =  {yl  +  v?) 12  for  a  single  blade  of  length  R  =77  mm,  rotating  in  the  plane  parallel  to 
the  ground  z  =  0  at  height  zo  =  R  at  frequency  85  Hz.  The  distance,  x,  is  measured  from  the 
axis  of  rotation.  The  circulation  of  the  tip  vortex  in  simulations  was  T  =  0.25m2/s. 

and  the  functionality  of  each  module  was  carefully  tested. 

One  of  the  major  developments  was  related  to  FMM  data  structures.  The  research  led  to  the 
development  of  a  novel  and  different  approach  to  the  computation  of  data  structures  that  avoids 
conventional  sorting  algorithms,  which  are  relatively  costly.  Instead,  algorithms  were  developed 
based  on  partial  sorting,  fixed-grid  sorting,  and  other  improvements  (see  list  of  publications).  Such 
improvements  enabled  accelerations  that  were  several  times  more  than  that  of  the  conventional 
algorithm  on  a  PC,  and  up  to  two  orders  of  magnitude  acceleration  on  the  GPU.  The  reason  why 
so  much  attention  was  placed  on  the  data  structures  and  this  part  of  the  FMM  algorithm  is  related 
to  the  fact  that  in  the  present  project  the  FMM  was  applied  to  dynamic  problems  where  particles 
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Figure  3.22:  Data  structure  construction  time  vs  number  of  particles  for  a  non-uniform 
(spherical-surface)  distribution  of  particles.  Computations  on  a  single  node  using  one  GPU. 
Each  sub-figure  corresponds  to  a  maximal  level  setting  for  the  octree. 

change  their  location  every  time  step,  and  efficient  construction  of  the  FMM  data  structures  be¬ 
comes  a  bottleneck,  especially  in  the  GPU  versions.  Figure  3.22  show  the  speedups  that  were 
obtained  for  the  construction  of  data  structures,  which  includes  both  hardware  and  algorithmic 
accelerations. 

One  of  the  achievements  of  the  present  research  is  a  full  FMM  on  a  single  GPU  suitable  for 
simulation  of  dynamics  problems.  The  very  first  implementation  of  the  FMM  on  the  GPU  did  not 
utilize  the  GPU  for  construction  of  data  structures,  as  it  was  found  to  be  relatively  inefficient  [7,  8]. 
This  algorithm  may  have  an  important  impact  on  a  number  of  fields  where  dynamic  FMM-based 
computations  are  important. 

Comparison  of  the  results  obtained  by  the  FMM  with  other  high-performance  techniques,  such 
as  the  FFT-based  algorithms  (pseudospectral,  DNS  simulations)  requires  solution  of  problems  with 
periodic  boundary  conditions,  which  within  the  known  approaches  requires  substantial  modifica¬ 
tion  of  the  core  FMM  code  and  data  structures.  Within  the  present  project  a  novel  method  to  obtain 
periodic  solutions  using  an  available  FMM  for  free-space  problems  is  developed  [13].  This  was 
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also  implemented  on  the  GPU.  Tests  demonstrate  a  good  performance.  The  algorithm  should  find 
applications  in  many  problems  in  CFD,  molecular  dynamics  and  crystallography. 

Heterogeneous  Algorithms  and  Architectures 

Despite  the  achievement  of  an  efficient  FMM  realization  on  a  single  GPU,  to  solve  larger 
problems  one  must  use  heterogeneous  architectures.  These  could  be  a  single  workstation  (with 
several  multicore  CPU  chips  and  a  few  GPUs),  or  a  distributed  network  of  such  inter-connected 
workstations. 

Careful  analysis  of  the  FMM  showed  that  some  parts  of  the  algorithm  cannot  be  efficiently 
mapped  on  GPUs  and  that  these  parts  have  comparable  or  better  performance  on  the  multicore 
CPU.  This  is  because  of  limited  local  memory  of  graphics  processors  (cache),  and  some  parts  of 
the  FMM  are  memory  bound.  However,  for  other  parts  the  GPU  can  be  very  efficient,  and  acceler¬ 
ate  them  by  a  factor  of  a  hundred  or  so.  Hence,  the  use  of  heterogeneous  CPU/GPU  architectures 
provides  a  natural  mapping/advantage  for  algorithms  with  different  parts  such  as  the  FMM.  We 
developed  this  concept  theoretically,  which  includes  parallel  execution  of  different  parts  of  the 
algorithm  on  the  CPU  and  GPU,  compacting  of  information  for  CPU/GPU  exchange,  data  struc¬ 
tures,  designing  for  different  levels  of  parallelism  (fine  and  coarse  grained),  and  developing  basic 
concepts  for  the  implementation  of  the  FMM  on  heterogeneous  architectures. 

A  paper  on  this  research  was  presented  at  the  Supercomputing’ll  conference  [4],  and  the  al¬ 
gorithm  has  since  been  subsequently  further  developed  and  extended  to  handle  vortex  elements 
[6,  12].  The  algorithm  showed  very  good  performance  and  scalability  and  enabled  computations 
of  test  problems  of  size  up  to  1  billion  particles  on  a  midsize  heterogeneous  cluster  (32  hetero¬ 
geneous  nodes  with  64  GPUs  on  the  Chimera  Cluster  at  the  University  of  Maryland  Institute  for 
Advanced  Computer  Studies).  While  different  performance  depending  on  the  problem  can  be 
achieved,  estimates  show  that  the  peak  performance  was  of  the  order  of  20  Teraflops,  but  if  the 
same  size  problem  should  be  computed  without  algorithmic  acceleration  (FMM),  a  2  Exaflop  per¬ 
formance  is  needed  to  achieve  the  observed  timing.  Figure  3.23  illustrates  timing  for  one  time  step 
achieved  by  this  algorithm  for  different  problem  sizes.  It  is  important  to  note  that  heterogeneous 
computing  is  a  trend  in  contemporary  high  performance  computing  because  it  utilizes  resources 
(multicore  PCs)  neglected  in  pure  GPU  computing,  provides  better  energy  efficiency,  substantial 
acceleration,  and  overall  economical  benefit. 
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Figure  3.23:  The  wall  clock  time  for  one  time  step  of  the  VEM  (vortex  blob  version)  acceler¬ 
ated  by  the  FMM  running  on  32  heterogeneous  nodes  of  the  UMIACS  Chimera  cluster.  The 
curves  correspond  to  the  case  of  a  scalar  potential  computation,  velocity  only,  and  velocity 
plus  stretching  term. 

Recent  Results 

During  this  past  year,  besides  continuation  tasks  as  development  of  user  friendly  interfaces  for 
the  developed  code  and  test  runs,  work  was  started  on  the  pilot  development  of  another  vortex 
element  method,  which  is  known  as  the  vortex-lattice  method  (VLM).  The  method  is  designed  for 
simulation  of  high  Reynolds  number  unsteady  flows  generated  by  lifting  surfaces,  and  has  been 
successfully  applied  to  solution  of  a  number  of  aerodynamics  problems.  A  feature  of  this  method 
is  that  it  computes  the  development  of  vortex  sheets,  as  opposed  to  vortex  filaments  and  blobs. 
In  principle,  such  an  approach  can  be  combined  with  other  vortex  methods,  such  as  Task  3.2. 
Accurate  description  of  the  flows  generated  by  rotorcraft  and  observed  in  experiments  requires 
modeling  and  computations  of  such  sheets;  see  Task  1.1  for  experimental  measurements. 

The  standard  implementation  of  the  VLM  is  even  more  limited  in  terms  of  the  size  of  the 
problem  than  the  vortex  particle  or  vortex  filament  methods  because  at  every  time  step  the  lifting 
surface  generates  additional  elements  of  the  sheet,  which  propagate  in  space  and  interact  with  solid 
objects.  After  several  time  steps  the  vorticity  distribution  becomes  very  complex,  the  number  of 
lattice  elements  very  large,  and  as  all  pairwise  interactions  need  be  taken  into  account,  making  the 
problem  computationally  very  challenging.  The  use  of  the  FMM  is  one  component  of  an  approach 
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Figure  3.24:  An  example  of  the  VLM  computations  accelerated  by  the  FMM.  Here  a  rect¬ 
angular  wing  rotates  and  creates  a  vortex  sheet  (ground  effect  is  taken  into  account,  which 
doubles  the  number  of  vertices  shown  on  the  upper  chart).  The  chart  on  the  bottom  shows  a 
slice  of  the  velocity  field  (xz  plane,  z  is  the  vertical  coordinate). 

to  reduce  the  problem  complexity.  Estimates  show  that  even  for  simple  surfaces  the  savings  at¬ 
tributed  to  the  FMM  can  be  several  orders  of  magnitude.  A  VLM  code  has  been  developed  and 
tested,  and  then  extended  to  use  the  FMM.  The  first  results  show  a  good  performance  and  poten¬ 
tial  for  its  use  in  the  future.  Figure  3.24  illustrates  a  preliminary  test  simulating  the  vortex  sheet 
generated  by  a  single  rotating  blade. 

The  VLM  code  then  was  used  as  an  example  to  create  user  friendly  Matlab  interface  for  hetero¬ 
geneous  computing  using  the  Parallel  computing  toolbox  (PCT).  Extensive  tests  of  the  VLM  using 
only  hardware  (GPU)  acceleration  and  both  algorithmic  and  hardware  acceleration  (heterogeneous 
FMM)  were  conducted.  While  successful  the  tests  also  showed  some  bottlenecks,  particularly, 
related  to  schemes  of  splitting  of  the  heterogeneous  algorithm,  which  should  be  addressed  in  the 
future. 
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Blade  Tip  Actuation  for  Tip  Vortex  Modifications  for  Brownout  Mitigation 
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Engineering,  the  University  of  Maryland,  and  The  U.S.  Naval  Academy 
Graduate  Student(s):  Julia  Elaine  Cline  and  Reza  Azizi  Ghasrehelal 
Contact  email(s):  dancila@uta.edu 

Background  and  Technical  Challenges 

Brownout  is  a  major  rotorcraft  operational  concern  in  military  rotorcraft  operations.  Brownout 
results  in  massive  particle  uplift  into  the  rotor  flow  and  can  lead  to  the  loss  of  visual  references  and 
situational  awareness  for  the  pilot  who  is  flying  the  aircraft.  Brownout  conditions  can  also  result 
in  excessive  erosion  of  rotor  blade,  engine  components  and  other  mechanisms. 

One  approach  to  mitigating  brownout  relies  upon  specific  flight  paths/maneuvers  in  ground 
proximity,  which  is  aimed  at  minimizing  the  quantity  of  particles  entrained  in  the  rotor  flow  and/or 
in  redirecting  those  already  entrained  away  from  the  field  of  view;  see  Task  3.3  for  details.  This 
type  of  mitigation  approach,  however,  may  restrict  the  flight  envelope  of  the  helicopter  and  cause 
operational  issues.  A  second  approach  relies  upon  electronic  ground  sensing  means  that  are  not 
impaired  by  brownout,  and  with  the  use  of  an  associated  cockpit  display  can  provide  the  crew  with 
sufficient  information  to  operate  the  aircraft  safely.  Both  of  these  approaches,  however,  do  not 
solve  the  excessive  abrasion  and  wear  caused  by  the  suspended  dust,  and  also  do  not  preserve  an 
unconstrained  operational  flight  envelope. 

Another  approach  to  brownout  mitigation,  and  the  approach  followed  in  this  task,  seeks  to 
significantly  reduce  particle  entrainment  into  the  rotor  flow.  This  approach  fundamentally  relies 
upon  the  recognition  that  while  brownout  is  a  complex,  two-phase  flow  phenomenon,  it  may  have 
a  flow  control  solution.  The  key  to  this  approach  is  an  understanding  of  the  fundamental  physical 
aspects  of  this  complex  flow.  To  this  end,  a  focused,  concerted  research  effort  in  this  direction  is 
being  undertaken  at  the  University  of  Maryland;  see  Tasks  1.1,  2.2  and  3.3.  The  complexity  of 
the  underlying  flow  physics  suggests  that,  ultimately,  a  combined  experimental  aerodynamics  and 
CFD  approach  is  the  most  profitable  avenue  of  research. 
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Technical  Approach 

Experimental  results  focused  on  brownout  studies  on  scaled  rotors  in  ground  effect  have  provided 
a  better  understanding  of  the  phenomenon  of  particle  entrainment  into  the  complex,  tip  vortex 
dominated  hovering  rotor  flow  in  ground  proximity;  see  Tasks  1.1  and  2.2.  The  results  show  that 
blade  tip  vortices  play  a  dominant  role  in  the  entrainment  of  particles  from  the  ground,  as  well  as  in 
their  subsequent  convection  into  the  rotor  flow.  The  entrainment  is  amplified  by  the  vortex  strength 
and  structure,  and  by  the  distance  between  the  vortex  and  the  ground.  This  observation  provides 
insight  into  possible  avenues  for  entrainment  mitigation,  i.e.,  increase  in  miss  distance  and/or  a 
reduction  in  vortex  strength  and/or  the  radial  diffusion  of  vorticity.  The  peak-to-peak  velocity  of  the 
vortex,  the  dominant  characteristic  for  particle  mobilization  and  entrainment,  is  not  determined  by 
vortex  strength  (circulation)  alone,  but  by  the  combination  of  vortex  strength  (circulation)  and  the 
vortex  core  structure  and  size  (vorticity  and  velocity  distribution).  Therefore,  a  strong  circulation 
vortex  with  a  larger  core  size  may  show  lower  particle  mobilization  and  entrainment  compared  to  a 
weaker  circulation  vortex  with  a  smaller  core  size.  This  behavior  suggests  that  it  should  be  possible 
to  reduce  the  entrainment  by  a  given  blade  tip  vortex,  with  a  given  circulation  level  (determined  by 
the  lift  developed  by  the  blade),  by  radially  diffusing  its  vorticity  and  increasing  its  core  size. 

The  objectives  of  this  research  are  to  investigate  the  validity  of  the  following  two  working  hy¬ 
potheses:  1.  Modulated  blade  tip  blowing — active  (piezoelectric)  or  passive  (resonant) — provides 
sufficient  flow  control  authority  to  modify  the  structure,  strength,  size,  and/or  trajectory  of  the 
blade  tip  vortex  such  that  with  proper  selection  of  input  the  particle  entrainment  capability  of  the 
blade  tip  vortex,  and  implicitly  that  of  the  rotor  flow  as  a  whole,  are  significantly  reduced,  resulting 
in  the  mitigation  of  brownout.  2.  Modulated  blade  tip  blowing — active  (piezoelectric)  or  passive 
(resonant) — provides  sufficient  flow  control  authority  to  excite  instabilities  of  the  tip  vortex  re¬ 
sulting  in  vortex  burst  and  consequently  spreading  vorticity  and  reducing  the  particle  entrainment 
capability  of  the  blade  tip  vortex,  and  implicitly  that  of  the  rotor  flow  as  a  whole,  resulting  in  the 
mitigation  of  brownout. 

Results 

The  focus  of  the  work  was  initially  upon  the  development  of  a  set  of  four  ducted  blade  configu¬ 
rations  representing  a  progression  from  the  baseline  slotted  blade  configurations  and  obtained  by 
migrating  the  blowing  orifices  to  the  upper  edge  and  by  progressively  varying  their  dihedral  orien¬ 
tation  in  increments.  These  blade  configurations  will  allow  a  comparative  investigation  of  steady 
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Figure  3.25:  A  modification  made  by  relocating  the  blowing  orifices  at  the  upper  edge  and 
oriented  at  45  degrees  dihedral.  Blade  tip  detail. 


Figure  3.26:  A  modification  made  by  increasing  the  dihedral  angle  to  60  degrees.  Blade  tip 
detail. 


Figure  3.27:  A  modification  made  by  reducing  the  dihedral  angle  to  30  degrees.  Blade  tip 
detail. 

blowing  effectiveness.  Four  blade  model  configurations  have  been  manufactured  that  implement 
steady  blowing  driven  by  centrifugal  pumping  and  with  blowing  nozzles  progressively  migrating 
from  the  baseline  configuration  (demonstrated  effective  previously  at  the  University  of  Maryland) 
towards  locations  and  orientations  that  may  provide  increased  effectiveness;  see  Figs.  3.25  to  3.27. 

Recent  work  on  this  task  has  been  performed  by  Milluzzo  and  Kuerbitz  at  the  U.S.  Naval 
Academy  under  the  Office  of  Naval  Research  (ONR),  award  number  N0001614WX30023,  exam¬ 
ining  the  fluid  dynamics  of  the  wake  produced  by  a  rotor  these  rotor  blades  as  it  interacted  with 
a  ground  plane.  Two-component  PIY  and  FV  was  used  to  evaluate  the  flow  field  generated  by  a 
rotor  operating  at  one  rotor  radius  above  the  ground.  The  goal  of  this  study  was  to  quantify  the 
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effect  centrifugal  pumping  blade  designs  had  on  the  flow  field  generated  by  the  rotor.  Five  different 
blades  were  tested;  a  rectangular  (non-slotted)  blade  was  used  as  a  baseline  case.  The  other  four 
blade  designs  had  an  inlet  slot  near  the  root  of  the  blade  and  a  internal  channel  that  connected  the 
inlet  slot  to  exit  slots  on  the  side  edge  of  the  blade.  The  exit  slots  were  designed  such  that  they 
exited  along  the  blade  chord  line  (i.e.,  0°  from  the  horizontal)  or  facing  the  upper  surface  of  the 
blade  at  30°,  45°,  and  60°  from  the  horizontal. 

Performance  measurements  were  also  taken  to  quantify  the  changes  in  performance  associated 
with  the  pumping  blade  designs.  Preliminary  results  have  indicated  that  at  lower  blade  collective 
pitch  angles  (i.e.,  lower  thrust)  the  rotor  with  the  baseline  blades  had  the  lowest  power  required. 
However,  as  the  rotor  thrust  was  increased  the  power  required  for  each  of  the  pumping  blade 
designs  began  to  approach  the  power  required  for  the  rotor  with  the  baseline  blades.  Furthermore, 
at  the  higher  values  of  thrust  the  blade  with  the  0°  exit  slots  required  less  power  than  for  the 
baseline  blades. 
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